Introduction The p63 protein encoded by the TP63 gene is a transcription factor of the p53 family and functions as a master regulator of ectodermal development. The key function of p63 during ectodermal development is underscored by phenotypic features in p63 knockout mice [1], [2] and in p63 knock-down zebrafish [3], [4]. The developmental abnormalities in animal models are reminiscent of those in p63-associated human disorders. Heterozygous mutations in p63 give rise to at least seven dominantly inherited clinical conditions with three major characteristics, ectrodactyly (also known as split hand/foot malformation, SHFM), orofacial clefting and ectodermal dysplasia with defects in skin, hair, teeth, nails and exocrine glands [5], [6]. There is a clear genotype-phenotype correlation in p63-associated disorders [7]. The most prominent of these disorders is the Ectrodactyly Ectodermal dysplasia and Cleft lip/palate syndrome (EEC, OMIM 604292) which combines all of the three phenotypic hallmarks and is almost invariably caused by missense mutations in the DNA binding domain of p63. Ankyloblepharon Ectodermal defects Cleft lip/palate syndrome (AEC, OMIM 106260) is caused by mutations in the SAM domain of the p63 that is involved in protein interaction. Nonetheless, mutations of p63 can explain only a minority of patients with only one of the three cardinal features, such as in patients with isolated SHFM (∼10%) and in patients with isolated cleft lip/palate (∼0.1%) [7]. There remains a large group of ectodermal dysplasia syndromes with phenotypes that resemble p63-associated syndromes [8]. The genetic basis of many of these clinically related conditions, referred to as the p63 phenotype network, is presently unknown. There is ample evidence that diseases clustering within such phenotype networks are caused by mutations in functionally related genes that constitute a gene network [9]–[11]. Elucidation of functional interactions among genes within the p63 gene network, their encoded proteins and regulatory elements which control expression of these genes will therefore provide new candidate genes for genetic disorders from the p63 phenotype network. Identifying target genes and cis-regulatory elements controlled by p63 is an important step in dissecting the p63 gene network. Previous studies have focused on transcriptional target genes of p63 identified through individual candidate gene approaches [12]–[14] or through genome-wide approaches [15]–[19]. However, the role of regulatory elements controlled by p63 in transcription has not yet been addressed so far. Split hand/split foot malformation (SHFM, OMIM 183600) is characterized by a deficiency of the central rays of the hands and feet, resulting in missing or malformed digits. SHFM may be isolated (non-syndromic) or be associated with other developmental anomalies (syndromic). Six distinct chromosomal loci for non-syndromic SHFM have been reported. Specific gene mutations have been identified in SHFM6 and SHFM4. SHFM6 (OMIM 225300, chromosome 12q13) is caused by a homozygous WNT10B mutation and it is the only autosomal recessive form of SHFM [20]. SHFM4 (OMIM 605289, chromosome 3q27) is caused by p63 mutations [21]. Chromosomal aberrations underlie three other types of isolated SHFM: 7q21 deletions and re-arrangements in SHFM1 (OMIM 183600) [22], 10q24 duplications encompassing the Dactylin gene (FBXW4) in SHFM3 (OMIM 600095) [23], and 2q31 deletions encompassing the HOXD gene cluster in SHFM5 (OMIM 606708) [24], [25]. In addition, linkage analysis has mapped SHFM2 (OMIM 313350) to chromosome Xq26 [26]. The SHFM1 locus on chromosome 7q21 has been delineated by various translocations, inversions, deletions and duplications [27]. The smallest region of overlapping deletions in SHFM1 patients [28] encompasses several genes: DYNC1I1, SLC25A13, DSS1, DLX5 and DLX6, of which only DLX5 and DLX6 have been shown to clearly play a role in early limb development. Dlx5 and Dlx6 are highly expressed in the apical ectodermal ridge (AER) of the developing limbs of mice [29]–[31] and in the fins of zebrafish [3], [4]. The AER is critical for limb outgrowth and patterning [32] and there is strong evidence that a failure to maintain the AER signaling is the main pathogenic mechanism in ectrodactyly [33]. The importance of the DLX5/DLX6 genes in limb development has been highlighted in mouse models. Dlx5 deficient mice do not show any limb defects [30]. However, an SHFM-like phenotype has been observed when both Dlx5 and Dlx6 were simultaneously deleted (Dlx5/Dlx6−/− ). The limb developmental phenotype in Dlx5/Dlx6−/− mice could be fully rescued by overexpression of Dlx5 in the AER [29], [34]. These observations suggest that the DLX5 and DLX6 genes cooperate in limb development by controlling a common developmental program. DLX5 and DLX6 are further expressed in the craniofacial prominence, the otic vesicle and in the brain [29]–[31], which correlates well with the hearing loss and mental retardation that are present in 30% of the SHFM1 patients [27]. While DLX5 and DLX6 are obvious candidate genes for SHFM1, mutations have not been found in either of the two genes. Here, we used a genome-wide DNA-binding profiling approach using Chromatin Immunoprecipitation (ChIP) followed by deep sequencing (ChIP-seq) in human primary keratinocytes to generate a catalogue of highly informative target genes and regulatory elements controlled by p63. One cis-regulatory element identified by DNA-binding profiling is located in the SHFM1 critical region and acts as an enhancer element for gene expression mediated by p63 during embryonic limb development. Our data indicate that loss of this element leads to SHFM1. This example illustrates that our catalogue of p63 binding sites can identify candidate genes and loci for the elucidation of disorders from the p63 phenotype network. Results Genome-wide p63 binding profile in human primary keratinocytes The most common isoform of p63, ΔNp63α, is highly expressed in the basal layer of the epidermis that consists mainly of keratinocytes. We therefore established human primary keratinocyte cultures (HKCs) from adult skin as our model system to elucidate p63 gene networks under physiological conditions. To identify target genes and regulatory elements controlled by p63, high-resolution global binding profiles of p63 were obtained from HKC cell lines established from two unrelated control individuals (wt1 and wt2) by ChIP-seq analysis using two antibodies recognizing different epitopes in p63 (4A4 and H129). Analysis of the sequenced reads using the peak recognition algorithm of Model-based Analysis of ChIP-Seq (MACS) [35] gave a highly significant overlap of 11,369 peaks from three profiles (P 1.5-fold compared to both backgrounds were determined as significant. We calculated the ROC AUC for all significant motifs and chose the best performing motif based on the ROC AUC (See Table S2 for the results). The PWM of this motif was combined with the p53scan algorithm to generate p63scan, using an optimal threshold, determined by the maximum f-measure as described previously [37]. The p63scan algorithm can be downloaded from http://www.ncmls.eu/bioinfo/p63scan/. To examine the correlation of motif score and peak height, all peaks were divided in quartiles according to peak height (the number of reads per peak). For each quartile the distribution of the motifs score as determined by p63scan is depicted as a boxplot. Motif over-representation analysis using TRANSFAC To detect putative transcription factor motifs reported in the TRANSFAC Professional database version 2009.3 [39], the MotifScanner program [91] was used. The search was performed on both strands using a 3rd-order Markov model calculated from the human promoter set of the Eukaryotic Promoter Database (EPD) as a background model. The parameter p (a prior probability of finding one instance of the motif in a sequence) was set to a value of 0.5. To identify motifs that are overrepresented in the p63 motif-less binding sites, the binomial test was used. The obtained P values were corrected for multiple testing (631 motifs for which sites were found in the p63-binding regions) using a Bonferroni correction. Quantitative PCR Quantitative PCR primers were designed using Primer3 (http://frodo.wi.mit.edu) [92], and qPCR reactions were performed in the 7500 Fast Real Time PCR System apparatus (Applied Biosystems) by using iQ SYBR Green Supermix (Biorad) according to the manufacturer's protocol. For qPCR of ChIP analysis, one primer set was used for each tested binding region (Table S8) and ChIP efficiency of certain binding sites was calculated using percentage of ChIPped DNA against input chromatin. Analysis of potential target genes associated with phenotypic defects using human and mouse disease bases Text mining-based [43] feature overrepresentation and gene to disease mapping were determined using the Online Mendelian Inheritance in Man (OMIM) disease database [93], [94]. Detailed information can be found in supplementary information. Human diseases associated with SHFM were taken from the Pictures Of Standard Syndromes and Undiagnosed Malformations (POSSUM) database [44], current as of August 2007, and mapped to genes through their OMIM IDs. Mouse SHFM-associated phenotypes and associated genes were taken from the Jackson Laboratory's Mouse Genome Database (http://www.informatics.jax.org/) [45]. Evolutionary conservation of binding sites To assess the evolutionary conservation of the 11,369 sites bound by p63, the PhastCons [54] conservation track from the UCSC Genome Browser was used to calculate PhastCons Conservation Score (PCCS). Conservation based on 44 vertebrate genomes was chosen because the p63 gene has 1-1 orthologs throughout the vertebrates [42]. The conservation for a region was calculated as the average conservation of each nucleotide therein. To analyse the correlation of PCCS and peak height, all peaks were divided in quartiles according to peak height (the number of reads per peak). For each quartile the distribution of the PhastCons Conservation Scores (PCCS) is depicted as a boxplot. Mapping deletion in a SHFM patient using ultra-high comparative genomic hybridization For detailed detection of chromosome 7 aberration, high resolution NimbleGen HG18 chromosome 7 specific 385K arrays were used (B3738001-00-01; Roche NimbleGen Systems, Madison, Wisconsin, USA). The 385K average probe distance was 365bp. DNA labeling, array hybridization, post-hybridization washes and scanning were performed according to the manufacturer's instructions (Roche NimbleGen). The acquired images were analyzed using NimbleScan V2.4 extraction software (Roche NimbleGen). For each spot on the array, the log2 Cy3/Cy5 ratio (relative intensity of the Cy3 labeled patient DNA vs. the Cy5 labeled male DNA reference pool of 5 healthy male individuals) was calculated using the segMNT algorithm, which also applied an automatic segment detection. A 50× averaging window was generated, resulting in 20kb segments for this array. Breakpoints were determined with SignalMap V1.9 software (Roche NimbleGen) and 20kb averaged log2 ratios were visualized in the UCSC genome browser. Constructs and transactivation assays The genomic regions of p63 binding site peaks were amplified by PCR with gateway cloning primers and cloned into a modified ccdB-containing pGL3-Enhancer Vector, or a ccdB-containing pGL3-Promoter Vector, or a ccdB-containing pGL3-Dlx5 Vector. The ccdB-containing pGL3-Dlx5 Vector was generated by amplification of mouse genomic DNA using primers described in Table S8 to obtain the mouse Dlx5 promoter to replace the SV40 promoter with BglII and HindIII sites in the ccdB-containing pGL3-Promoter Vector. Point mutations were introduced into p63-binding motifs of SHFM1-BS1 to generate mutant p63 binding sites, where the essential cytosine and guanine bases were mutated to adenosine. The ΔNp63α wild-type (pcDNA-mM_ΔNp63α) expression plasmid has been described previously [85]. Point mutations were introduced into this plasmid to generate R204W, R279H, R304W, K194E and L517F mutations. Transfection and luciferase assays were described previously [85]. All cloning and mutagenesis primers are described in Table S8. Transgenic reporter analyses in zebrafish and mice Human genomic fragments containing the SHFM1-BS1, -BS2 and -BS3 were amplified with primers described in Table S8. The PCR fragments were subcloned in PCR8/GW/TOPO vector and then transferred, through recombination using Gateway technology, to the ZED destination vector for zebrafish transgenesis [95]. This vector contains the Xenopus Cardiac actin promoter driving DsRed as a positive control for transgenesis. To generate the zebrafish transgenic embryos, we used Tol2 transposon/transposase method [96] with minor modifications. Volume of 2–5nl of mixture containing 25ng/ul of transposase mRNA, 20ng/ul of phenol/chloroform purified ZED constructs and 0.05% phenol red was injected in the cell of 1 cell stage embryos. Three or more independent stable transgenic lines were generated for each construct. For the generation of transgenic mice, the genomic fragments with and without point mutations in p63 consensus motif were transferred into a vector containing the human minimal beta-globin promoter, lacZ and a SV40 polyadenylation signal. Constructs were linearized and the vector backbone removed prior to microinjection into the pronucleus of one-cell mouse embryos. F0 embryos of 9.5–13 dpc stages were harvested and stained for lacZ activity. Morpholino injections and in situ hybridizations in zebrafish embryos Once cell stage embryos were injected with 3ng of ΔNp63 MO II (TCCACAGGCTCCAGGATTCTTACCC) as described previously [4]. Injected embryos were raised at 28°C in standard E3 medium and fixed at 48 hours post fertilization in 4% paraformaldehyde overnight at 4°C. In situ hybridizations were carried out as described [97]. As a control, we injected a similar amount of a MO directed against the Xenopus tropicalis olig2 gene that shows no match in the zebrafish genome [98]. Chromosome Conformation Capture (3C) assay Chromosome Conformation Capture (3C) assay was performed as referred in Hagege et al., 2007 [56]. Limbs of E10- and E15-stage mouse embryos were dissected and processed to get single cells preparations. Ten million isolated cells were first fixated with 2% formaldehyde, and then cells were lysed and nuclei were digested with HindIII endonuclease (Roche). After that, DNA was ligated with T4 DNA ligase (Promega) in low concentration conditions to favour intramolecular ligations. A set of locus specific primers close to a HindIII site (Table S8) was designed with Primer3 v. 0.4.0 [92]. These primers were used to make semi-quantitative PCRs to measure the relative enrichment in each ligation product. The primer near to the BS1 enhancer was taken as the fixed primer, and the different interactions were tested using primers close to the promoters of DLX5 and DLX6 genes. For each interaction two negative control primers were designed about 30 kb upstream and downstream the promoter specific primer. PCR products were run in agarose gels and measured using a Typhoon scanner. Product values were related to a control composed of two BACs that encompass our region of interest. Supporting Information Figure S1 ChIP-qPCR analysis of p63 binding in human primary keratinocytes using two different p63 antibodies 4A4 (pan-p63) and H129 (α-specific). Specific binding of p63 to the tested binding sites was observed, including to binding sites at p21WAF/CIP19 and DST which served as positive controls, but not to the negative controls myoglobin exon 2 (myo) and a no-gene region (chr11). (0.16 MB PDF) Click here for additional data file. Figure S2 Correlation of motif scores to peak heights. All peaks were divided in quartiles according to peak height (the number of reads per peak). For each quartile the distribution of the motifs score as determined by p63scan is depicted as a boxplot. (0.16 MB PDF) Click here for additional data file. Figure S3 Genomic distribution of p63 binding sites with and without p63 consensus binding motifs. (0.13 MB PDF) Click here for additional data file. Figure S4 Correlation of PhastCons Conservation Scores (PCCS) to peak heights and p63 binding motifs of the p63 binding sites. (A) The percentage of p63 binding sites (y-axis) is plotted against decreasing cut-off values of the PhastCons Conservation Score (PCCS) (x-axis) for two groups of peaks: those with a p63 motif and those without a p63 motif, as determined by p63scan. (B) All peaks were divided in quartiles according to peak height (the number of reads per peak). For each quartile the distribution of the PhastCons Conservation Scores (PCCS) is depicted as a boxplot. (0.08 MB PDF) Click here for additional data file. Figure S5 Previously reported chr.7 deletions involved in SHFM1 (hg18). (A) The grey track CGH array CNV data with deletion in the SHFM1 patient is compared to the minimum critical regions for SHFM1 based on the literature. Pale green tracks represent deleted intervals in patients with SHFM1 who have cytogenetic deletions (Del Porto et al., 1983; Tajara et al., 1989; Morey et al., 1990; Roberts et al., 1991; Nunes et al., 1994; McElveen et al., 1995; Montgomery et al., 2000). Dark green tracks represent deleted intervals in patients with SHFM1 where mapping has been done with STS markers (Marinoni et al., 1995; Crackower et al., 1996; Fukushima et al., 2003; Wieland et al., 2004). Purple tracks represent summed mapping of deletions combined from many patients with SHFM1 (Scherer et al., 1994; Tackels-Horne et al., 2001). The brown track represent a microdeletion at the break point of a chromosome inversion in a patient with SHFM1 (Brown et al., 2010). The red bar represents three p63 binding sites. (B) A zoomed-in view of the region including SHFM1-BS1, -BS2 and -BS3 and DLX5/6. A translocation in SHFM1 that disconnects SHFM1-BS1 with DLX5/6 is depicted with a black arrow (Saitu et al., 2009). (0.24 MB PDF) Click here for additional data file. Figure S6 Genomic qPCR analysis of deletion of the p63 binding site SHFM1-BS1 in the SHFM1 patient. Genomic qPCR was performed to confirm the deletion in the SHFM1 patient revealed by CGH analysis. Copy number was calculated against an internal control primer set in the CFTR gene. (0.18 MB PDF) Click here for additional data file. Figure S7 Raw data of 3C experiments in mouse limbs at E10 and E15. (0.47 MB PDF) Click here for additional data file. Figure S8 Evolutionary conservation of the p63 binding site SHFM1-BS1 in vertebrates. SHFM1-BS1 was examined for its conservation in vertebrates using USCS genome browser. The consensus motif of p63 is highlighted in the red box. (0.19 MB PDF) Click here for additional data file. Table S1 Validation by ChIP-qPCR of detected binding sites in ChIP-seq analysis. (0.07 MB DOC) Click here for additional data file. Table S2 Motif analysis with a de novo motif prediction pipeline. (0.02 MB XLS) Click here for additional data file. Table S3 The motifs significantly overrepresented in the p63 motif-less binding sites. (0.03 MB DOC) Click here for additional data file. Table S4 The most significant GO annotation terms of potential target genes of p63. (0.05 MB DOC) Click here for additional data file. Table S5 OMIM IDs associated with genes containing p63 binding sites. (0.14 MB XLS) Click here for additional data file. Table S6 Affected features in p63 potential target gene-associated diseases (PTG-associated diseases). (0.37 MB XLS) Click here for additional data file. Table S7 Genes potentially involved in SHFM mouse models or human SHFM phenotypes. (0.04 MB DOC) Click here for additional data file. Table S8 Primers for expression, ChIP analyses, and for cloning. (0.03 MB XLS) Click here for additional data file.
