Genome-Wide DNA Methylation Analysis During Palatal Fusion Reveals the Potential Mechanism of Enhancer Methylation Regulating Epithelial Mesenchyme Transformation

Abstract

Epithelial mesenchyme transformation (EMT) of the medial edge epithelium (MEE) is the crucial process during palatal fusion. This work is aimed to elucidate the enhancer regulatory mechanism by genome-wide DNA methylation analysis of EMT during palatal fusion. Over 800 million clean reads, 325 million enzyme reads, and 234 million mapping reads were generated. The mapping rate was 68.85–74.32%, which included differentially methylated 17299 CCGG sites and 2363 CCWGG sites (p < 0.05, log₂FC >1). Methylated sites in intron and intergenic regions were more compared to other regions of all DNA elements. GO and KEGG analysis indicated that differential methylation sites related to embryonic palatal fusion genes (HDAC4, TCF7L2, and PDGFRB) at the enhancer were located on CCWGG region of non-CpG islands. In addition, the results showed that the enhancer for HDAC4 was hypermethylated, whereas the enhancers for TCF7L2 and PDGFRB were hypomethylated. The methylation status of enhancer regions of HDAC4, PDGFRB, and TCF7L2, involved in the regulation of the EMT during palatal fusion, may enlighten the development of novel epigenetic biomarkers in the treatment of cleft palate.

Introduction

Cleft palate (CP), which may account for 30% of cleft lip and/or palate in humans, is a congenital malformation caused by the failure of palate shelves to oppose at the midline and fuse (Stuppia et al., 2011; Rahimov et al., 2012). It is universally acknowledged that the epithelial mesenchyme transformation (EMT) of the medial edge epithelium (MEE) is the most crucial process in embryonic palatal fusion. However, regulated mechanisms at the molecular level, especially the regulation of enhancer methylation during embryonic palatal development, are poorly understood. In the mouse, palatal shelves appose at the midline and palatal fusion occurs at an embryonic gestation day of 14.5. It was previously reported that, during the process of EMT, the unbalance of embryonic palate mesenchyme (MEPM) cell proliferation and apoptosis might result in CP (Rice, 2005; Thiery and Sleeman, 2006; Nawshad, 2008).

DNA methylation that occurs at the C5 position of cytosine within CpG and non-CpG (CpA, CpC, and CpT) in the genome is regarded as a critical component of epigenetic regulation (Barres et al., 2009; Law and Jacobsen, 2010). Moreover, DNA methylation plays a crucial role in many biological processes, such as embryogenesis, cellular differentiation, X-chromosome inactivation, genomic imprinting, and transposon silencing (Jaenisch and Bird, 2003; Levin and Moran, 2011; Peters, 2014), and has attracted attention because of its wide range of action, reversibility, and genetic characteristics. Methylation of CpG islands in gene promoter regions correlates with transcriptional silencing (Caiafa and Zampieri, 2005), by either inhibition of transcription factor binding (Geiman and Robertson, 2002) or interaction of methyl-CpG binding proteins associated with transcriptional repression (Deaton and Bird, 2011). Recent studies showed that methylation of CpG islands within the enhancer regions greatly influenced malignant transformation and progression (Aran and Hellman, 2013; Aran et al., 2013), which suggests that enhancers of CP-related genes (HDAC4, TCF7L2, and PDGFRB) in embryonic mice may be a new clue to epigenetic markers. However, patterns of methylation of CpG islands within the enhancer regions of genes underlying palate development, and their contribution to CP are still unclear. Moreover, identifying genes whose enhancers undergo differential methylation during palate development would be a significant step toward predicting and deciphering the pathogenesis of CP. Our work aims to elucidate the differentially methylated enhancers and regulatory mechanisms of EMT during palatal fusion by genome-wide DNA methylation analysis.

Although most of the epigenetic markers are involved in methylation of target promoter regions, more and more evidence suggests that methylation in enhancer regions correlates better with target gene expression than promoter methylation (Aran et al., 2013). Enhancers are 20-400 bp DNA sequences, located in CpG islands or non-CpG islands, which combine with tissue-specific transcription factors and regulate gene transcription at distant loci through chromosome looping (Petrascheck et al., 2005). Enhancers of their activity can be highly tissue specific, and the tight interplay between active enhancers and binding of their specific transcription factors (Feldmann et al., 2013; Kieffer-Kwon et al., 2013) are intricately dependent on enhancer methylation (Saluz et al., 1986). Hypomethylated enhancers allow transcriptional activation of target genes. Conversely, hypermethylated enhancers usually result in the silencing of target genes (Shlyueva et al., 2014). Nevertheless, details on the enhancer methylation patterns of EMT during palatal fusion are very limited. To explore the mechanisms of CP development during palatal fusion, we sequenced the embryonic mouse palatal tissues (E14.5), using MethylRAD combined with high-throughput sequencing (HiSeq X Ten), and generated the first genome-wide DNA methylation profile of embryonic mouse palatal tissues. Simultaneously, we screened the enhancer methylation levels of genes associated with the formation of CP, and subjected the genes with differentially methylated sites to Gene Ontology (GO) and KEGG pathway enrichment analysis.

HDAC4 has a central role in the formation of the skeleton (Vega et al., 2004), and controls the development of the palatal skeleton (Haberland et al., 2009). Hypermethylated HDAC4 can inhibit gene expression and transcription factor binding to target DNA during palatal fusion, and may lead to the formation of CP. TCF7L2 overexpression can inhibit EMT and palatal fusion. PDGFRA and PDGFRB interaction during palatogenesis in both zebra fish and mice had affected proper condensation of maxillary neural crest cells during palate formation (McCarthy et al., 2016). Hypomethylated PDGFRB can promote PDGFRB overexpression and inhibition of EMT during palatal fusion, contributing CP formation. These findings illustrated that HDAC4, TCF7L2, and PDGFRB are the candidate genes associated with palate development, and any abnormal expression of HDAC4, TCF7L2, and PDGFRB may lead to CP during palatal shelve fusion.

Subjects and Methods

Ethics, animals, and treatment

This study was approved by The Animal Ethics Committee of Shantou University Medical College (Guangdong, China). C57BL/6J mice were purchased from the Beijing Vital River Laboratory Animal Technology Co. Ltd. (Beijing, China), and were housed in the animal experimental center of Shantou University Medical College. Females were mated with fertile males overnight with similar weight and ages, and the embryonic day 0.5 (E0.5) was designated at 8 AM of the next day when a vaginal plug was found. Pregnant females at E10.5 were randomly divided into two groups. The case group was administered with at-RA (Sigma-Aldrich, MO) that was dissolved in corn oil, using gavage at 70 mg/kg, whereas the control group was given the equivalent volume of corn oil according to their weight. C57BL/6J mice at E14.5 were killed and palatal shelves were dissected using microbiological instruments. The palatal shelves were immediately frozen in liquid nitrogen and then kept at −80°C for further research.

MethylRAD library preparation and sequencing

Genomic DNA was extracted from embryonic mouse palatal shelves using the conventional cetyltrimethylammonium (cetrimonium) bromide (CTAB) method, and the genomic library was prepared using MethylRAD (in presence of an activator; D = A or G or T; S = C or G) (Cohen-Karni et al., 2011; Wang et al., 2015). Base quality values were calculated using a Phred quality score (Q sanger = −10log₁₀ p). Input sequencing data before operation and computing were called raw data or raw reads. After operation, the data, including adapter reads and low-quality sequences, were removed from raw reads as clean reads. The clean reads were then subjected to Pair-End sequencing on a HiSeq X Ten platform by the Shanghai Oebiotech Co. Ltd (Shanghai, China). To improve the accuracy in the follow-up analysis, Pear software (v0.9.6) (Zhang et al., 2014) was used to refilter pair-end sequencing paired clean reads according to the following terms: (i) remove low quality reads (more than 15% of base mass lower than 30) and (ii) remove sequences containing too many N bases (N containing more than 8%). After mapping of the paired-end sequencing reads, quality control was conducted. Finally, the MethylRAD sequence data were obtained while the labels without target enzyme sites were filtered and deleted based on each sample's location in the library. Base distribution and quality distribution maps of clean reads were plotted (Fig. 1). [Reference genome: ftp://ftp.ensembl.org/pub/release-84/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.toplevel.fa.gz.Annotation information: ftp://ftp.ensembl.org/pub/release-84/gff3/mus_musculus/Mus_musculus. GRCm38.84.chr.gff3.gz). (Fig. 2: schematic overview of the procedure for MethylRAD library preparation (Wang et al., 2015)].

FIG. 1.

(A) Base percentage composition (the abscissa represents the pair-end reads position, reflecting the proportion of A/C/G/T/N bases in each position, and the ordinate represents percent of A/C/G/T/N bases). (B) Distribution of base qualities (the abscissa represents the pair-end reads position, and the ordinate represents the base quality values).

FIG. 2.

Schematic overview of the procedure for MethylRAD library preparation. Genomic DNA is digested with the restriction enzyme FspEI, producing 32-bp fragments, including four base 3’overhangs. Adaptors with compatible overhangs (NNNN) are ligated to each end of these fragments. Tag density can be adjusted using adaptors with selective overhangs (e.g., NNNG). The constructs are amplified and purified by gel extraction. Sample-specific barcodes are incorporated in each construct by PCR, and the products pooled for sequencing.

Distribution of DNA methylation sites

According to the genome sequence, a label containing CCGG/CCWGG sites was proposed as the reference sequence, and the individual enzyme reads were mapped to the reference sequence of the CCGG/CCWGG sites using the SOAP program (version 2.21, parameter: -M4-v2-r0) (Li et al., 2009). A reliable methylation site was determined by finding the sequencing depth of no less than 3 sites. Considering that there was a repetition region in the reference sequence, reads in the repetition region were simultaneously compared to a plurality of positions in the reference sequence. It was difficult to judge the actual position of the methylation site.

Using MethylRAD technology combined with high-throughput sequencing, the degree of methylation at a site could be determined by observing the sequencing depth of the methylation site. The density distribution map of methylation sites depicts the distribution density of CCGG/CCWGG sites and the actual detected methylation site density on each chromosome. Detecting the distribution of sites on the chromosome group using MethylRAD technology was more uniform. Moreover, by using the window sliding method to calculate the CCGG/CCWGG number on each window on the chromosome, the site distribution over the entire chromosome was described (Supplementary Fig. S1; Supplementary Data are available online at www.liebertpub.com/dna). The sequencing depths of the methylation sites (CCGG/CCWGG) in each sample were plotted using a box plot (the number of methylation sites in each sample had a depth higher than 3) (Fig. 3).

FIG. 3.

Box plots showing the sequencing depths of the methylation sites (CCGG and CCWGG) in each sample (case (b1, b2, and b3) and control (B1, B2, and B3) embryonic palatal tissues). In the box chart, the box part is the main body of the box-shaped chart, and the middle of the black horizontal line is the median of the data; the upper and lower sides of the box is a quarter of the data that are greater than the upper quartile (Q3), and a quarter of the data are less than the lower quartile (Q1). The interval between Q1 and Q3 is called the interquartile range (IQR). The longitudinal lines in the upper and lower sides of the box are tentacle lines. The upper cutoff line of the tentacle line is “Q3 + 1.5 * IQR” and the lower one is “Q1-1.5 * IQR”.

The UTR region was calculated using snpEff software (version: 4.3p) (Cingolani et al., 2012), and was counted using bed tools software (v2.25.0) (Quinlan and Hall, 2010) according to the annotation document and the distribution of methylation sites in the different gene elements in each sample. The gene elements included the 3′UTR, 5′UTR, TSS2000, exon region, intron region, and intergenic region. The distribution of methylation sites on different functional elements of the genome is shown in Figure 4.

FIG. 4.

Distribution in different components of the genome. (A) CCGG methylation sites. (B) CCWGG methylation sites. The y-axis shows the number of methylation sites. The x-axis shows the different components of the genome.

Relative quantification of DNA methylation levels

According to the consistency of the equal-length sequence amplification efficiency, the sequencing depth of the methylation label can reflect the methylation level of the site (CCGG/CCWGG). The RPM (reads per million) is the unit of the quantitative value of methylation level at the site, which means the site methylation level quantitative value equal to site coverage reads number/library high quality reads multiplied by 1,000,000. The DNA methylation level in the gene region was calculated by the reads that were aligned on a distinct locus in the genome. We selected upstream and downstream of 2 kb sections of the gene body, TSS and TTS, and summarized the DNA methylation level of the distribution trend of sequencing reads in embryonic palatal shelve (Supplementary Fig. S2).

Comparison of methylation levels between groups

The change in methylation level was assessed based on the sequencing depth information of each site in the relative quantitative results of methylation, using R package edge R (Robinson et al., 2010). The differential p-value (p < 0.05) and fold change (log₂FC >1) between different sites were given according to the sequencing data (Supplementary Data). For assessing the methylation level of differential methylation sites between case and control for the three biological replicates, we implemented cluster analysis to further reveal the changes in CCGG/CCWGG methylation levels between the groups. The cluster heat map is shown (p-value of differential methylation sites selected by the top 10,000) in Figure 5.

FIG. 5.

Hierarchical cluster analysis heat-map of differential methylation sites between groups. (A). The heatmap of differential CCGG methylation sites, (B) the heat-map of differential CCWGG methylation sites.

According to the positional information of the differentially methylated sites relative to the chromosome, the distribution Circos of differential methylation sites on chromosomes was drawn and the results are shown in the Figure 6.

FIG. 6.

Distribution of differential methylation sites on chromosomes ((A) CCGG differential methylation site; (B) CCWGG differential methylation site.) The sites of upregulated and downregulated methylation were expressed in dark gray and light gray, respectively.

The pie map distribution of differential methylation sites on different functional components was drawn according to the positional information of the differentially methylated sites relative to the gene, and the results are shown in the Figure 7. The distribution of differential methylation sites on genes of different functional components showed that the CCGG/CCWGG is mainly in the intron and the intergenic regions.

FIG. 7.

Distribution of differential methylation sites on genes of different functional components, (A) CCGG differential methylation site; (B) CCWGG differential methylation site.

GO and KEGG enrichment analysis

Differentially methylated sites related to genes involved in the formation of CP were screened for GO enrichment analysis. The numbers of genes included in each GO category were counted and the statistical significance of gene enrichment in each GO category was calculated using the hypergeometric distribution test. GO entries with the number of corresponding genes greater than 2 in three categories were screened and the GO enrichment analysis results (Top 30) are shown in Figure 8.

FIG. 8.

GO enrichment analysis of differentially methylated genes. (A) CCGG GO enrichment analysis; (B) CCWGG GO enrichment analysis, including cellular component, molecular function, and biological process. GO, gene ontology.

We also performed KEGG enrichment pathway analysis between the groups (mapping to annotate genes) and screened differentially methylated sites for KEGG enrichment analysis, and identified differentially methylated genes participating in cell signaling pathways (Fig. 9).

FIG. 9.

KEGG enrichment of top 20 bubble diagram of differential methylation-related genes. (A) CCGG KEGG analysis, (B) CCWGG KEGG analysis. The x-axis is the enrichment score, the bigger the bubble, the more genes they contain. The color of bubbles varied from 2e–04, 4e–04, 6e–04 to 8e–04, and the concentration indicated by p-value is as large as the enrichment score. The x-axis comprises the cell signal pathways.

Identifying potential enhancers

The BLAST Software (McGinnis and Madden, 2004) (www.ncbi.nlm.nih.gov/) was used to analyze the enriched enhancer region combined with the information of annotation gene among the hypermethylated and hypomethylated genes in the case and control (Supplementary Fig. S3). The difference and similarity of enhancer in two groups were analyzed. The common enhancers were regarded as the final enhancer of the target genes.

Statistical analysis

For the paired case and control of embryonic mouse palatal tissue, an unsupervised hierarchical clustering analysis was used to identify distinct subgroups based on the differentially methylated sites. Relative quantification of DNA methylation levels was performed using the normalized read depth (RPM) for each site. The correlation between samples of methylation levels was assessed using Pearson's correlation coefficient. The level of methylation between groups was assessed using R package edge R (Robinson et al., 2010). A p-value <0.05 and log₂FC >1 were considered statistically significant.

Results

Identification of the distribution of methylated sites

We generated a total of over 800 million clean reads from six embryonic mouse palatal shelf samples. Of the total reads, 325 million enzyme reads and 234 million mapping reads were generated. The mapped ratio between cases and controls was 68.85–74.32% (Table 1), which included differentially methylated 17299 CCGG sites and 2363 CCWGG sites that were mapped to the reference sequence. The average sequencing depth of the methylation sites (CCGG/CCWGG) in each sample between the case and control embryonic palatal shelf is shown in Table 2.

Table 1.

Sample Sequencing Data Quantity Versus Ratio (The Ratio of Enzyme Reads to the Unique Position of the Reference Sequence Is Shown; B1, B2, and B3 Are Cases, b1, b2, and b3 Are Controls)

Sample	Clean data	Enzyme reads	Mapping reads	Mapped ratio
B1	131,029,920	52,959,691	38,391,986	72.49%
B2	131,029,920	56,224,181	41,472,605	73.76%
B3	131,029,920	54,365,634	38,595,473	70.99%
b1	142,120,402	59,873,387	38,595,473	74.32%
b2	142,120,402	58,784,358	42,436,660	72.19%
b3	131,029,920	42,987,516	29,595,826	68.85%

Table 2.

The Average Sequencing Depth of the Methylation Site (CCGG and CCWGG) in Each Sample (Case (B1, B2, and B3) and Control (b1, b2, and b3) Embryonic Palatal Tissues)

Sample	CCGG sites number	Mean depth	CCWGG sites number	Mean depth
B1	752,150	38.6	57,545	9.07
B2	815,022	37.12	50,898	8.33
B3	770,026	37.65	55,469	8.32
b1	823,508	38.29	78,578	8.15
b2	790,142	40.32	69,122	8.65
b3	767,888	29.16	57,252	8.36

Analysis of the methylation peaks among the chromosomal distribution of the CCGG/CCWGG site revealed that the number of methylation peaks among cases was greater compared with controls in most chromosomes (Supplementary Fig. S1). The distribution of methylation sites on different gene elements in each sample indicates that the onset of cases during palatal shelf fusion might be associated with global methylation: (i) methylated sites (CCGG/CCWGG) in intron and intergenic regions were more than other regions of all DNA elements; and (ii) among all genomic functional components, the methylation sites in the cases were significantly higher than those in the controls. Except for 5′UTR in CCGG, the methylation sites in the case group were significantly lower than those in the controls (Tables 3 and 4).

Table 3.

The RPM of Different Methylation Sites in HDAC4, PDGFRB, and TCF7L2 That Are Related to Embryonic Palatal Fusion Genes That Regulate EMT in Each Sample

Gene	Position	B1RPM	B2RPM	B3RPM	b1RPM	b2RPM	b3RPM
HDAC4	92052147	5.42	5	9.02	3.25	1.82	0
PDGFRB	61071357	2.17	1.25	1.13	4.06	7.27	5.44
TCF7L2	55852561	1.08	2.5	0	3.25	7.27	5.44

Cases (B1, B2, and B3) and Controls (b1, b2, and b3) embryonic palatal tissues.

EMT, epithelial mesenchyme transformation; RPM, reads per million.

Table 4.

Position and Sequence of Different Methylation Sites and Stat-End of Enhancer Region in HDAC4, PDGFRB, and TCF7L2

Gene	Position	Sequence	Stat-End of enhancer region	Log2FC	p-value
HDAC4	92052147	GAAGCTCAGACCGCCAGGAGGGTGCACTCAA	92051600-92053400	1.7	0.011
PDGFRB	61071357	TGGGGTGGCATGGCCTGGTCCCCAGGGTGCC	61070801-61071800	−1.8	0.007
TCF7L2	55852561	GAGGTGGGTGGGCCCTGGCCTGTTTACTCCA	55850001-55852600	−2.1	0.005

Identification of methylation level

The results of DNA methylation level showed that (i) the DNA methylation level of the distribution trend of sequencing reads in cases were higher than controls; (ii) the DNA methylation level in gene body, TSS and TTS of CCGG/CCWGG, was dramatically increased in the 2 kb region upstream, except for gene body of CCWGG; (iii) the peak of DNA methylation level of CCGG was located on 300 bp in gene body of 1 kb in TSS and −1 kb in TTS, whereas the peak of DNA methylation level of CCWGG was located on start site in gene body of 200 bp in TSS and 600 bp in TTS; and (iv) the DNA methylation level in gene body and TTS of CCGG/CCWGG was dramatically dropped in the transcription termination site, whereas TSS of CCGG/CCWGG was dramatically increased and dropped to the lowest point before the TSS. However, MethylRAD technology is dependent on high sequences, such that only regions of relatively high CCGG/CCWGG abundance are obtained (Supplementary Fig. S2). Apparently, most of the genes with spaced methylation sites may not be detected by MethylRAD, even if few methylation sites present are actually methylated.

Hierarchical cluster analysis results indicate that the methylation levels of differential methylation sites among cases are higher than those of controls. Hypomethylated CCGG/CCWGG in cases and controls are clustered near the bottom, whereas hypermethylated CCGG/CCWGG in cases and controls are clustered near the top. Interestingly, we found that there were distinctive methylation patterns in cases and controls, and showed obvious interindividual and intraindividual differences in methylation profiles among groups (Fig. 5). The differential methylation sites on chromosomes of the brown track represented upregulated genes, whereas the blue track represented downregulated genes (Fig. 6).

Analysis of distribution of methylation sites in CCGG/CCWGG among genes of the different functional components showed the genome coverage in exon (9.54%), intergenic (39.4%), upstream (6.68%), 3′ UTR (2.07%), 5′ UTR (0.51%), and intron (42.82%) in CCGG; and exon (11.3%), intergenic (35.26%), upstream (7.1%), 3′ UTR (2.85%), 5′ UTR (0.67%), and intron (42.82%) in CCWGG (Fig. 7).

Screening for potential enhancer of gene related to CP

To evaluate differentially methylated genes that are related to CP during palatal shelf fusion, the enhancer within genes (HDAC4, PDGFRB, and TCF7L2) among the sequences of differentially methylated sites in the cases and controls was screened based on enhancer annotation and the method of sequence alignment (Tables 5 and 6). We thoroughly scrutinized the list of significantly differentially methylated loci for sites located within genes (HDAC4, PDGFRB, and TCF7L2) that are related to embryonic palatal fusion and responsible for regulating EMT. The position of differentially methylated site of HDAC4, PDGFRB, and TCF7L2 is located within CCWGG, non-CpG islands of intron region (Supplementary Fig. S3) (www.ncbi.nlm.nih.gov/). The differentially methylated site of HDAC4 enhancer was hypermethylated, whereas PDGFRB and TCF7L2 were hypomethylated.

Table 5.

Number of Distribution of Methylation CCGG Sites in the Different Gene Elements

	B1	B2	B3	b1	b2	b3
Exon	90310	97981	92752	98042	95125	92224
Intergenic	334952	363396	342433	367520	351298	341827
Intron	378475	409860	387778	414481	398146	386672
Upstream	58900	65032	60581	64611	62211	59871
Utr3prime	19997	21481	20549	21792	21050	20440
Utr5prime	3363	4339	3522	3767	3620	3356

Table 6.

Number of Distribution of Methylation CCWGG Sites in the Different Gene Elements

	B1	B2	B3	b1	b2	b3
Exon	8407	7650	7912	10674	9563	7069
Intergenic	22570	19785	22198	31477	27567	24998
Intron	31104	27637	29738	42535	37388	29288
Upstream	5722	5050	5568	7784	6776	5481
Utr3prime	1845	1673	1746	2370	2124	1590
Utr5prime	461	396	489	637	561	435

GO and KEGG enrichment analysis for HDAC4, PDGFRB, and TCF7L2

GO enrichment analysis was performed to identify differentially methylated genes related to HDAC4, PDGFRB, and TCF7L2 and the results of enrichment (p < 0.001) were observed for GO biological process and molecular function, which are shown in Table 7. Up on KEEG pathway enrichment analysis, we further identified differentially methylated genes related to HDAC4, PDGFRB, and TCF7L2, and the results of enrichment are shown in Table 8.

Table 7.

Gene Ontology Enrichment Analysis of HDAC4, PDGFRB, and TCF7L2 in Cases and Controls of Palate Shelve

Gene	Biological process	Molecular function
HDAC4	Regulation of striated muscle cell differentiation; regulation of skeletal muscle fiber development	Protein kinase binding; RNA polymerase III transcription factor binding;
PDGFRB	Cell migration involved in coronary angiogenesis; metanephric glomerular mesangium development; protein autophosphorylation; cell migration; signal transduction	ATP binding; protein kinase binding
TCF7L2	Catenin import into nucleus; regulation of hormone metabolic process	Protein kinase binding; nuclear hormone receptor binding

Table 8.

KEEG Pathway Enrichment Analysis of HDAC4, PDGFRB, and TCF7L2 (p < 0.001) in Cases and Controls of Palate Shelve

Gene	KEEG pathway	p value	Enrichment score
HDAC4	Human papillomavirus infection	0.00090	1.73
PDGFRB	microRNAs in cancer	0.00002	2.55
	Focal adhesion	0.00007	2.17
	Gap junction	0.00035	2.63
	Glioma	0.00070	2.66
	Human papillomavirus infection	0.00090	1.73
	Regulation of actin cytoskeleton	0.00100	1.89
	Ras signaling pathway	0.00130	1.83
TCF7L2	Human papillomavirus infection	0.00092	1.73
TCF7L2	Pathways in cancer	0.00397	1.53

Discussion

In this study, we performed a genome-wide DNA methylation analysis, after generating a genome-wide map of DNA methylation during palatal fusion that regulates EMT. Moreover, we further identified the genes with several differentially methylated sites and analyzed the pathways by comparison with case-matched control. The most striking finding in this analysis is the distinct and characteristic methylation pattern displayed in the enhancer region of intron, as reflected by the number of differentially methylated non-CpG sites and corresponding genes.

EMT plays a crucial part in the embryogenesis of palatal fusion. Particularly, in mesoderm formation, neural crest cell delamination, fibrosis, and wound healing, EMT is indispensable (Kalluri and Weinberg, 2009; Thiery et al., 2009). In the palate formation, the apposing shelves adhere to each other by adhesion junctions, then attached by desmosomes in the MEE cell. During EMT, the MEE seems disappearing, leads to the confluence of connective tissue across the palate. CP has a multifactorial etiology that often includes failure of adherence of opposing individual palatal shelves and/or EMT of the MEE (Nawshad, 2008). Failure of any of these evolutionarily late, finely orchestrated processes occurs relatively often and results in CP (Rice, 2005).

A considerable amount of data has established that enhancers are not merely a collection of TF binding sites but also have the capacity to drive transcription independence of their target promoters. This feature of enhancers suggests that they serve more regulatory functions than previously appreciated (Kim and Shiekhattar, 2015). Andersson et al. (2014, 2015) showed that enhancers are enriched in core promoter elements, but are CpG poor, whereas promoters are generally CpG rich. Beside the CpG content, enhancers and promoters have broad similarities and overlapping functional properties, and have been considered to form a single class of regulatory element (Andersson et al., 2015). Therefore, the non-CpG enhancer of intron region participates in the regulation of gene expression. Gene expression is precisely regulated by the combination of promoters and gene-distal enhancer regulatory regions (Vernimmen and Bickmore, 2015).

Aberrant DNA methylation of the susceptible genes, especially regulatory region can result in gene silencing; therefore, enhancers play pivotal roles in regulating gene expression during embryonic palatal fusion. Methylation-induced gene silencing can occur either by impeding the binding of transcription factors and/or the methyl group serving as a substrate for methyl-CpG binding proteins (Bogdanović and Veenstra, 2009). The latter may, in turn, recruit chromatin remodeling proteins, which affect various histone modifications to establish a repressive chromatin environment that results in dramatically altered levels of gene expression (Deaton and Bird, 2011).

Tissue-specific DNA methylation patterns are precisely programmed during embryogenesis (Kafri et al., 1992). Failure to establish correct methylation patterns may lead to embryonic lethality (Li et al., 1992) or result in developmental craniofacial malformations (Ohgane et al., 2001), including CP (Loenarz et al., 2010). All-trans retinoic acid (at-RA), an endogenous metabolite of vitamin A, is required for normal pattern formation during embryogenesis (Ackermans et al., 2011). Conversely, abnormally high concentrations in both experimental animals and humans result in fetal malformations, including CP (Cuervo et al., 2002; Ackermans et al., 2011). Campbell et al.'s (Campbell et al., 2004) research shows the maximum effect for limb defects, and total malformed fetuses were seen at 60 mg/kg and higher. The incidence of CP is increased over the entire range of administered doses, reaching a maximum of 74% (100 mg/kg). They have concluded that a single 100 mg/kg of at-RA to mice bearing E11 embryos has been shown to result in a 90–100% malformation frequency in mouse embryos. In our previous experimental work, the incidence of CP in the at-RA 70 mg/kg and 100 mg/kg treatment groups was 100%, and the incidence of CP in the at-RA30 mg/kg group was 40%. A rate of CP as assessed by our work appears consistent with the results of all earlier studies, that doses of 70 mg/kg or higher induce CP (Qin et al., 2014). Therefore, this study used 70 mg/Kg doses of at-RA to induce the CP model of embryo mice.

HDAC4 [histone deacetylase 4] encoding histone deacetylase 4 is a class II histone deacetylase that binds to other HDACs and myocyte enhancing factor-2 (Mef2) to inhibit transcription factor binding to target DNA (Haberland et al., 2009). Recent studies have suggested that HDAC4 has a central role in the formation of the skeleton (Vega et al., 2004) and controls the development of the palatal skeleton (Haberland et al., 2009). In this study, the differentially methylated non-CpG site of HDAC4, located in the enhancer region, was hypermethylated, which can inhibit gene expression and transcription factor binding to target DNA during palatal fusion, and may lead to the formation of CP.

TCf7L2 [transcription factor 7-like 2 (T cell specific, HMG-box)] plays a critical role in the Wnt signaling pathway (Sun et al., 2010). A previous study suggested that TCF7L2 is a member of the TCF lymphocyte enhancer factor (LEF) protein family (MacDonald et al., 2009), and the bipartite transcription factor β-catenin/TCF-LEF serves as an effector of cAMP-dependent protein kinase A (PKA) signaling. The TCF7L2-enhanced signaling mediates the physiological effects of peptide hormones, including glucagon-like peptide-1 (GLP-1), which utilizes cAMP as a second messenger (Jin et al., 2008; Liu and Habener, 2008). Yet, TCF4 (officially known as TCf7L2) is a candidate gene associated with CP (Essaoui et al., 2013). In this study, the differentially methylated non-CpG site of TCF7L2, located in the enhancer region, was hypomethylated in CP mice, and likely promotes TCF7L2 overexpression and negative regulation of the canonical Wnt signaling pathway to inhibit EMT and palatal fusion.

PDGFRB [platelet-derived growth factor receptor, beta polypeptide] encodes a receptor tyrosine kinase. Recent studies have suggested that PDGFRB can regulate epithelial-to-mesenchymal transition of epicardial cells (Smith et al., 2011). In addition, while loss of PI3K signaling through PDGFRB results in normal craniofacial development, the combined loss of PI3K signaling through PDGF receptors may result in craniofacial defects, as severe as a PDGFRA single mutant (Klinghoffer et al., 2002). Moreover, PDGFRA and PDGFRB are shown to interact during palatogenesis in both zebra fish and mice, and affects proper condensation of maxillary neural crest cells during palate formation (McCarthy et al., 2016). In our study, the differentially methylated non-CpG site of PDGFRB, located on the enhancer region, was hypomethylated, which can promote PDGFRB overexpression, and negatively regulates the PI3K-Akt signaling pathway to inhibit EMT during palatal fusion, leading to the formation of CP.

Although we obtained several enhancers showing differential methylation that regulates EMT during palatal fusion, we do not exclude that other factors may drive CP formation, nor are we suggesting a direct causality between enhancer methylation and CP; some limitations of our study must be taken into account. First, the sample size was relatively small. Second, palatal shelves were obtained directly from embryonic mouse and may be mixed with other tissues. Third, enhancers exhibit highly dynamic methylation patterns during embryo development, whereas our study is only the process of embryonic palatal fusion and lacks of the whole development of the embryonic palatal shelf. Fourth, due to DNA methylation is maintained by the activity of methyltransferase (Burgers et al., 2002), we cannot exclude the effect of at-RA and its endogenous small molecule metabolites on methyltransferase, which further affect the bias of experimental results. To evaluate the enhancer DNA methylation in complex tissues, laser capture microdissection must be employed, and the DNA methylation site and level of gene must be thoroughly verified. Therefore, DNA methylation site and gene-level verification need further evaluation and could be focused as an intending research direction.

Conclusions

Enhancer methylation leading to transcriptional silencing is a major reason for the inactivation of genes. Our results demonstrate that EMT during palatal fusion is dependent on the regulation of enhancer methylation status. The methylation status of enhancer regions of HDAC4, PDGFRB, and TCF7L2 well provides a foundation for future studies. The enhancer of HDAC4, PDGFRB, and TCF7L2 may serve as a future epigenetic target for development of novel therapeutics for CP.

Footnotes

Acknowledgment

This work was supported by The National Natural Science Foundation of China Project grant (number 81001284), Guangdong Province Natural Science Foundation Project grant (number 2015A030313431), and Shanghai Oebiotech Co. Ltd (Shanghai, China).

Disclosure Statement

The authors have declared that no competing interests exist.

References

Ackermans

, Zhou

, Carels

, Wagener

, and Von den Hoff

(2011). Vitamin A and clefting: putative biological mechanisms. Nutr Rev, 69, 613–624.

Andersson

, Gebhard

, Miguel-Escalada

, Hoof

, Bornholdt

, Boyd

, et al. (2014). An atlas of active enhancers across human cell types and tissues. Nature, 507, 455–461.

Andersson

, Sandelin

, and Danko

(2015). A unified architecture of transcriptional regulatory elements. Trends Genet, 31, 426–433.

Aran

, and Hellman

(2013). DNA methylation of transcriptional enhancers and cancer predisposition. Cell, 154, 11–13.

Aran

, Sabato

, and Hellman

(2013). DNA methylation of distal regulatory sites characterizes dysregulation of cancer genes. Genome Biol, 14, R21.

Barres

, Osler

M.E.

, Yan

, Rune

, Fritz

, Caidahl

, et al. (2009). Non-CpG methylation of the PGC-1alpha promoter through DNMT3B controls mitochondrial density. Cell Metab, 10, 189–198.

Bogdanović

, and Veenstra

(2009). DNA methylation and methyl-CpG binding proteins: developmental requirements and function. Chromosoma, 118, 549–565.

Burgers

, Fuks

, and Kouzarides

(2002). DNA methyltransferases get connected to chromatin. Trends Genet, 18, 275–277.

Caiafa

, and Zampieri

(2005). DNA methylation and chromatin structure: the puzzling CpG islands. J Cell Biochem, 94, 257–265.

10.

Campbell

, Smith

, Fisher

, and Warren

(2004). Dose-response for retinoic acid-induced forelimb malformations and cleft palate: a comparison of computerized image analysis and visual inspection. Birth Defects Res B Dev Reprod Toxicol, 71, 289–295.

11.

Cingolani

, Platts

, Wang le

, Coon

, Nguyen

, Wang

, et al. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly, 6, 80–92.

12.

Cohen-Karni

, Xu

, Apone

, Fomenkov

, Sun

, Davis

P.J.

, et al. (2011). The MspJI family of modification-dependent restriction endonucleases for epigenetic studies. Proc Natl Acad Sci U S A, 108, 11040–11045.

13.

Cuervo

, Valencia

, Chandraratna

, and Covarrubias

(2002). Programmed cell death is required for palate shelf fusion and is regulated by retinoic acid. Dev Biol, 245, 145–156.

14.

Deaton

A.M.

, and Bird

(2011). CpG islands and the regulation of transcription. Genes Dev, 25, 1010–1022.

15.

Essaoui

, Nizon

, Beaujard

M.P.

, Carrier

, Tantau

, de Blois

M.C.

, et al. (2013). Monozygotic twins discordant for 18q21.2qter deletion detected by array CGH in amniotic fluid. Eur J Med Genet, 56, 502–505.

16.

Feldmann

, Ivanek

, Murr

, Gaidatzis

, Burger

, and Schubeler

(2013). Transcription factor occupancy can mediate active turnover of DNA methylation at regulatory regions. PLoS Genet, 9, e1003994.

17.

Geiman

T.M.

, and Robertson

K.D.

(2002). Chromatin remodeling, histone modifications, and DNA methylation-how does it all fit together?. J Cell Biochem, 87, 117–125.

18.

Haberland

, Montgomery

R.L.

, and Olson

E.N.

(2009). The many roles of histone deacetylases in development and physiology: implications for disease and therapy. Nat Rev Genet, 10, 32–42.

19.

Jaenisch

, and Bird

(2003). Epigenetic regulation of gene expression: how the genome integrates intrinsic and environmental signals. Nat Genet, 33 Suppl, 245–254.

20.

Jin

, George Fantus

, and Sun

(2008). Wnt and beyond Wnt: multiple mechanisms control the transcriptional property of beta-catenin. Cell Signal, 20, 1697–1704.

21.

Kafri

, Ariel

, Brandeis

, Shemer

, Urven

, McCarrey

, et al. (1992). Developmental pattern of gene-specific DNA methylation in the mouse embryo and germ line. Genes Dev, 6, 705–714.

22.

Kalluri

, and Weinberg

(2009). The basics of epithelial-mesenchymal transition. J Clin Invest, 119, 1420–1428.

23.

Kieffer-Kwon

K.R.

, Tang

, Mathe

, Qian

, Sung

M.H.

, Li

, et al. (2013). Interactome maps of mouse gene regulatory domains reveal basic principles of transcriptional regulation. Cell, 155, 1507–1520.

24.

Kim

, and Shiekhattar

(2015). Architectural and functional commonalities between enhancers and promoters. Cell, 162, 948–959.

25.

Klinghoffer

R.A.

, Hamilton

T.G.

, Hoch

, and Soriano

(2002). An allelic series at the PDGFalphaR locus indicates unequal contributions of distinct signaling pathways during development. Dev Cell, 2, 103–113.

26.

Law

J.A.

, and Jacobsen

S.E.

(2010). Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet, 11, 204–220.

27.

Levin

H.L.

, and Moran

J.V.

(2011). Dynamic interactions between transposable elements and their hosts. Nat Rev Genet, 12, 615–627.

28.

, Bestor

, and Jaenisch

(1992). Targeted mutation of the DNA methyltransferase gene results in embryonic lethality. Cell, 69, 915–926.

29.

, Yu

, Li

, Lam

T.W.

, Yiu

S.M.

, Kristiansen

, et al. (2009). SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics (Oxford, England), 25, 1966–1967.

30.

Liu

, and Habener

J.F.

(2008). Glucagon-like peptide-1 activation of TCF7 L2-dependent Wnt signaling enhances pancreatic beta cell proliferation. J Biol Chem, 283, 8723–8735.

31.

Loenarz

, Ge

, Coleman

, Rose

, Cooper

, Klose

, et al. (2010). PHF8, a gene associated with cleft lip/palate and mental retardation, encodes for an Nepsilon-dimethyl lysine demethylase. Hum Mol Genet, 19, 217–222.

32.

MacDonald

B.T.

, Tamai

, and He

(2009). Wnt/beta-catenin signaling: components, mechanisms, and diseases. Dev Cell, 17, 9–26.

33.

McCarthy

, Liu

J.S.

, Richarte

A.M.

, Eskiocak

, Lovely

C.B.

, Tallquist

M.D.

, et al. (2016), Pdgfra and Pdgfrb genetically interact during craniofacial development. Dev Dyn, 245, 641–652.

34.

McGinnis

, and Madden

T.L.

(2004). BLAST: at the core of a powerful and diverse set of sequence analysis tools. Nucleic Acids Res, 32, W20–W25.

35.

Nawshad

(2008). Palatal seam disintegration: to die or not to die?. that is no longer the question. Dev Dyn, 237, 2643–2656.

36.

Ohgane

, Wakayama

, Kogo

, Senda

, Hattori

, Tanaka

, et al. (2001). DNA methylation variation in cloned mice. Genesis, 30, 45–50.

37.

Peters

(2014). The role of genomic imprinting in biology and disease: an expanding view. Nat Rev Genet, 15, 517–530.

38.

Petrascheck

, Escher

, Mahmoudi

, Verrijzer

C.P.

, Schaffner

, Barberis

(2005). DNA looping induced by a transcriptional enhancer in vivo. Nucleic Acids Res, 33, 3743–3750.

39.

Qin

, Shen

, Peng

, Wu

, Hu

, Zhang

, et al. (2014). Metabolic characterization of all-trans-retinoic acid (ATRA)-induced craniofacial development of murine embryos using in vivo proton magnetic resonance spectroscopy. PLoS One, 9, e96010.

40.

Quinlan

A.R.

, and Hall

I.M.

(2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics (Oxford, England), 26, 841–842.

41.

Rahimov

, Jugessur

, and Murray

J.C.

(2012). Genetics of nonsyndromic orofacial clefts. Cleft Palate Craniofac J, 49, 73–91.

42.

Rice

D.P.

(2005). Craniofacial anomalies: from development to molecular pathogenesis. Curr Mol Med, 5, 699–722.

43.

Robinson

M.D.

, McCarthy

D.J.

, and Smyth

G.K.

(2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (Oxford, England), 26, 139–140.

44.

Saluz

H.P.

, Jiricny

, and Jost

J.P.

(1986). Genomic sequencing reveals a positive correlation between the kinetics of strand-specific DNA demethylation of the overlapping estradiol/glucocorticoid-receptor binding sites and the rate of avian vitellogenin mRNA synthesis. Proc Natl Acad Sci U S A, 83, 7167–7171.

45.

Shlyueva

, Stampfel

, and Stark

(2014). Transcriptional enhancers: from properties to genome-wide predictions. Nat Rev Genet, 15, 272–286.

46.

Smith

C.L.

, Baek

S.T.

, Sung

C.Y.

, and Tallquist

M.D.

(2011). Epicardial-derived cell epithelial-to-mesenchymal transition and fate specification require PDGF receptor signaling. Circ Res, 108, e15–e26.

47.

Stuppia

, Capogreco

, Marzo

, La Rovere

, Antonucci

, Gatta

, et al. (2011). Genetics of syndromic and nonsyndromic cleft lip and palate. J Craniofac Surg, 22, 1722–1726.

48.

Sun

, Wang

, and Jin

(2010). Insulin alters the expression of components of the Wnt signaling pathway including TCF-4 in the intestinal cells. Biochim Biophys Acta, 1800, 344–351.

49.

Thiery

, Acloque

, Huang

, and Nieto

(2009). Epithelial-mesenchymal transitions in development and disease. Cell, 139, 871–890.

50.

Thiery

J.P.

, and Sleeman

J.P.

(2006). Complex networks orchestrate epithelial-mesenchymal transitions. Nat Rev Mol Cell Biol, 7, 131–142.

51.

Vega

R.B.

, Matsuda

, Oh

, Barbosa

A.C.

, Yang

, Meadows

, et al. (2004). Histone deacetylase 4 controls chondrocyte hypertrophy during skeletogenesis. Cell, 119, 555–566.

52.

Vernimmen

, and Bickmore

(2015). The hierarchy of transcriptional activation: from enhancer to promoter. Trends Genet, 31, 696–708.

53.

Wang

, Lv

, Zhang

, Dou

, Sun

, Li

, et al. (2015). MethylRAD: a simple and scalable method for genome-wide DNA methylation profiling using methylation-dependent restriction enzymes. Open Biol, 5, pii: 150130.

54.

Zhang

, Kobert

, Flouri

, and Stamatakis

(2014). PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics (Oxford, England), 30, 614–620.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.25 MB

0.28 MB

0.17 MB

0.97 MB