Genome-wide DNA methylation profile of developing deciduous tooth germ in miniature pigs

Background DNA methylation is an important epigenetic modification critical to the regulation of gene expression during development. To date, little is known about the role of DNA methylation in tooth development in large animal models. Thus, we carried out a comparative genomic analysis of genome-wide DNA methylation profiles in E50 and E60 tooth germ from miniature pigs using methylated DNA immunoprecipitation-sequencing (MeDIP-seq). Results We observed different DNA methylation patterns during the different developmental stages of pig tooth germ. A total of 2469 differentially methylated genes were identified. Functional analysis identified several signaling pathways and 104 genes that may be potential key regulators of pig tooth development from E50 to E60. Conclusions The present study provided a comprehensive analysis of the global DNA methylation pattern of tooth germ in miniature pigs and identified candidate genes that potentially regulate tooth development from E50 to E60. Electronic supplementary material The online version of this article (doi:10.1186/s12864-016-2485-9) contains supplementary material, which is available to authorized users.


Background
Tooth development is characterized by a sequential and reciprocal series of inductive signals transmitted between the epithelium and mesenchyme [1,2]. A number of transcription factors, signaling molecules, growth factor receptors, and extracellular matrix molecules have been identified as participating in the regulation of odontogenesis patterning and differentiation processes [1]. Currently, almost all of our understanding of molecular mechanisms controlling tooth formation and mineralization has come from studies in mice [3][4][5], whose teeth are significantly different from those of humans in regards to both morphology and number, with only one set of dentition and an absence of canines and premolars [6]. Recently, our group used miniature pigs, which resemble humans in anatomy, physiology, pathophysiology, and development, as an animal model to study the complicated mechanism of tooth development. We identified the characteristic patterns of the spatiotemporal morphogenesis of successional teeth in miniature pigs [7] and mapped the mRNA and microRNA expression profiles [8,9]. Our results identified differential gene expression patterns in different developmental stages and a spatio-temporal pattern of down-regulation during tooth development [8], suggesting that tooth formation is the result of tight control by a sequence of molecular networks that act at particular places and times. Despite numerous studies leading to the discovery of genetic mechanisms underlying the expression of specific genes during different stages of tooth development [3], recent studies suggest that epigenetic mechanisms also participate in tooth development [10][11][12].
DNA methylation is one of the best-studied epigenetic modifications and regulates a variety of processes, including embryonic development, cellular differentiation, tissue-specific gene expression, genomic imprinting, X chromosome in activation, and chromosome stability [13][14][15]. In vertebrates, DNA methylation occurs predominantly at cytosine residues within CpG dinucleotides [16]. CpG-rich regions are referred to as CpG islands (CGIs), and are located in nearly 40 % of the promoters of mammalian genes [17]. CGIs usually have global unmethylated patterns [18]. Unmethylated CGIs in the promoter regions are normally associated with gene expression, whereas methylated CGIs usually result in gene silencing [13,19]. Although DNA methylation has been widely accepted as a key mechanism of transcriptional regulation and a critical factor in the development of various organs [20][21][22][23], little is known about the normal developmental changes in DNA methylation during odontogenesis.
In the present study, we characterized the genome-wide temporal dynamics of DNA methylation in developing deciduous molars in miniature pigs using methylated DNA immunoprecipitation combined with high-throughput sequencing (MeDIP-seq) to further investigate the role of DNA methylation in gene expression during tooth development. We revealed the landscape of the methylome and identified the potential role of DNA methylation in gene expression in developing tooth germ. The results support utilizing the pig as a model organism for tooth development research.

Results
Global mapping of DNA methylation in the tooth germ of miniature pigs We mapped the global DNA methylation status of tooth germ collected from E50 and E60 miniature pig embryos. After the removal of low-quality reads from raw MeDIPseq data, an average of 4.6 Gb of clean reads were obtained per sample. Total reads were mapped to the reference genome, and mapping rates ranged from 84.30 to 85.55 %, of which 59.95-62.85 % were uniquely mapped to specific regions in the pig genome (Additional file 1: Table S1). Only uniquely mapped reads were used in further analysis.
Uniquely mapped reads were detected in all chromosomes (GGA1-18 and chromosome X) (Additional file 2: Figure S1). The genome coverage of the CG, CHG, and CHH sites negatively correlated with read depth; most regions had low sequencing depth, and a small number of regions had highs equencing depth (Additional file 3: Figure S2). The distribution of reads indifferent CG density regions was also analyzed, demonstrating that low CG density regions covered more uniquely mapped reads than any other region (Additional file 4: Figure S3). Different genomic regions exhibited different methylation patterns. The majority of reads were present in the intron regions, followed by CGIs (Fig. 1a). Within the gene body, a depletion of or decrease in reads occurred at the transcription start site (TSS) in both E50 and E60 tooth germ. In contrast, a gradual increase in reads occurred in the intragenic region (Fig.1b).
Methylation peaks, referred to as methyl-cytosine-enriched regions, are important parameters for the identification of global DNA methylation status [24]. Using MACS1.4.0., we obtained 135,043-179,112 peaks, covering 6.97-7.94 % of the reference genome (Additional file 1: Table S2). We further analyzed the distribution of peaks in different components of the genome. A major portion of the peaks were present in the CDS region, followed by 2 kb downstream of the transcription termination site and 2 kb upstream of the TSS, whereas the 5'-UTR, intron, and 3'-UTRhad fewer peaks (Fig. 2). Here, 2 kb upstream of TSS was considered the proximal promoter.

MeDIP-seq data validation by bisulfite sequencing
To validate the MeDIP-seq data, three regions with relatively high methylation were selected at random for bisulfite sequencing. The results obtained for the three gene regions were in accordance with the MeDIP-seq results. Moreover, great consistency was found among individuals in the same group (Fig. 3, Additional file 5: Figure S4).

Analysis of functional categories of differentially methylated genes
To identify the biological functions associated with differentially methylated genes at E60 versus E50, we performed Gene Ontology (GO) analysis using the DAVID program. GO assignments revealed that the up-methylated genes were mainly involved in 301 categories (p < 0.05, Additional file 7). The top 20 GO categories was showed in Fig. 5a including signal transduction, proteolysis, extracellular matrix organization, small molecule metabolic processes, cell-matrix adhesion, axon guidance, and intracellular signal transduction (p < 3E-06). The downmethylated genes strongly correlated with 180 categories (p < 0.05, Additional file 8). The top 20 included transcription (DNA-dependent), negative regulation of apoptosis, regulation of transcription (DNA-dependent), signal transduction, positive regulation of transcription (DNAdependent), transport, and negative regulation of chromatin silencing at rDNA (p < 0.001, Fig. 5b).

Pathway and path-net analysis of differentially methylated genes
To determine the significant pathways involved in differential methylation, we used the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database to predict putative functions. The up-methylated genes were significantly enriched in 79 predicted pathways (p < 0.05, Additional file 9), the most significant of which was regulation of the actin cytoskeleton, inositol phosphate metabolism, proteoglycans in cancer, phosphatidylinositol signaling, and focal adhesion (p < 0.0005, Fig. 5c). The down-methylated genes were tightly related to 18 pathways, including gastric acid secretion, endocytosis, the wnt signaling pathway, and cell cycle (p < 0.05, Fig. 5d, Additional file 10).
To further understand the interactions among the significant pathways involving the up-and down-regulation of differentially methylated genes, an interaction net was built using path-net analysis. As shown in Fig. 6a, some differentially methylated genes were identified as being involved in key pathways during tooth morphogenesis, including the calcium signaling pathway, apoptosis, focal adhesion, pathways in cancer, regulation of the actin cytoskeleton, and ErbB signaling.

Signal-net analysis
Gene signal transduction networks (Signal-net) were built to explore inter-gene signaling between the differentially methylated genes. A total of 104 genes were screened as potential core regulators of tooth development from E50 to E60 (Additional file 11). As shown in Fig. 6b, PLCG2, PLA2G4A, PLD1, PIK3CB, PIK3CD, JAK2, PPAP2C, MAPK10 had higher betweenness centrality, degree, indegree, and outdegree, indicating that they may play crucial roles in odontogenesis. GSK3β, ERBB4, and IGF1R, which are involved in the wnt signaling pathway, ErbB signaling, and insulin-like growth factor (IGF) signaling, respectively, were also identified as core regulators that interact with other molecules in the signaling network.

Discussion
This study provides a comprehensive analysis of genomewide DNA methylation profiles in the developing tooth germ of miniature pigs. The results revealed analogous DNA methylation patterns to those of other species, from mammals to plants [24][25][26][27][28][29][30][31]. In the tooth germ of miniature pigs, the CGIs and promoters (around TSSs) remain hypomethylated, whereas the methylation levels in gene body regions were relatively high. In mammals, most CpG sites are methylated, and highly methylated sequences are usually found in the gene body, repeats, satellites, and non-repetitive intergenic regions [32]. However, CGIs in approximately 60 % of the promoters of human genes and TSSs are normally hypomethylated [33]. In general, the DNA methylation levels of CGIs and TSSs are associated with gene expression [34,35]. Unmethylated CGIs in the promoters are essential for tissue-specific expression of corresponding genes in early embryos and later somatic cells [36]. The methylation patterns revealed in this study were consistent with those reported in previous studies.
The comparison of DNA methylation profiles between E60 and E50 tooth germ showed the different DNA methylation patterns in the different developmental stages of pig tooth germ. Pathway analysis of differentially methylated genes indicated that many tooth morphogenesisrelated pathways were involved. As expected, several significant pathways were obtained, including those involved in the actin cytoskeleton, focal adhesion, wnt signaling, and calcium signaling. Among these pathways, the wnt pathway is one of the most important during tooth morphogenesis. From E50 to E60, pig tooth germ progressed from the bell stage to secretory stage. At the bell stage, the enamel-and dentine-forming cells differentiate. Initially, the mesenchymal cells differentiate into odontoblasts that secrete the dentine matrix. Subsequently, the adjacent epithelial cells differentiate into ameloblasts and secrete enamel matrix. The odontoblasts and ameloblasts then control the phases of mineralization of enamel and dentine. At the secretory stage, tooth development is characterized by the transformation of soft gel-like extracellular matrices into mineralized structures of the skeleton [37]. The actin cytoskeleton, focal adhesion, and calcium signaling pathway play key roles in cytoskeletal organization [38]; the pathway analysis implied that they likely contribute to tooth biomineralization. The pathway profiles involving differentially methylated genes provided insight into understanding the maturation mechanisms of tooth germ in miniature pigs.
Signal-net analysis screened several key genes that may contribute to the transformation of soft tissue into Fig. 2 Distribution of DNA methylation peaks in different genomic regions. Peaks were defined as regions with sequencing tags more than 10 and p value < 10 −5 with MACS. A major portion of peaks were present in the CDS region, followed by 2 kb downstream from the transcription termination site and2 kb upstream from the TSS, whereas the5'-UTR, intron, and 3'-UTRhad fewer peaks. The x-axis shows the different genomic regions. The y-axis is the number of peaks hard tissue, including PLCG2, PIK3CB, PIK3CD, GSK3β, MAPK10, ERBB4, and IGF1R. The phospholipase C gamma-2 (PLCG2) gene is important for intracellular signal transduction pathways because it encodes an enzyme that plays a crucial role in the generation of second messengers following the hydrolysis of phosphatidylinositol 4, 5-bisphosphate [39]. GSK3β, a negative regulator of the wnt canonical pathway, plays a crucial role in tooth morphogenesis, and inhibition of GSK3β could delay the differentiation of ameloblasts and odontoblasts [40]. PIK3CB and PIK3CD belong to the phosphoinositide 3-kinase (PI3K) family. PI3K and MAPK family members are involved in the regulation of many cellular processes, such as proliferation, migration, survival, and apoptosis, as well as tooth development [41][42][43]. ERBB4 and IGFIR were also identified to be associated with odontogenesis and ameloblast differentiation [44][45][46]. The results indicate novel epigenetic mechanisms regulating tooth mineralization.

Conclusions
In summary, this study provided a comprehensive analysis of DNA methylation profiles in developing tooth germ from miniature pigs during biomineralization and identified 104 differentially methylated genes that may be potential core regulators of tooth development from E50 to E60. The genes and pathways screened suggest strong candidates for in-depth studies of the epigenetic mechanisms underlying tooth development in miniature pigs.

Animals and tissue collection
The Wuzhishan miniature pigs used in this study were purchased from Kexing Laboratory Animal Company of Beijing, China. This study was conducted in accordance with the recommendations of the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China). All the animal experiments were approved by the Animal Care and Use Committees of Capital Medicine University under the permit (Approval Number CMU-B20100106). The last deciduous mandibular molar germs at embryonic day 50 (E50) and E60 which corresponded to early and late bell stages were used in this study. Pregnant miniature pigs were anesthetized and sacrificed as previously described [8]. The tooth germs were isolated under a microscope, frozen immediately in liquid nitrogen, and stored separately at−80°C until DNA extraction.

MeDIP sequencing
For each developmental stage, the tooth germs were isolated separately from three randomly selected embryos with a specific sex (male) as biological replicates. MeDIP DNA libraries were constructed following the protocol as described previously [47]. In brief, genomic DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Five microgram DNA was sonicated to fragments ranging from 100 to 500 bp. Subsequently, DNA

Bioinformatic analysis
Raw sequencing data were first processed to filter out low-quality reads that containing more than 5 'N's or over 50 % of the sequence with low quality value (Phred score < 5). The MeDIP-seq data were aligned to the UCSC pig reference genome (Sscrofa9.2, http://hgdownload.cse.ucsc.edu/goldenPath/susScr2/chromosomes/.), allowing up to two mismatches using SOAP2 (Version 2.21) [48]. For each sample, we analyzed the distributions of MeDIP-Seq reads on the pig genome, in CGIs, gene region and different CG density regions, as well as the genome coverage of the CG, CHG, and CHH (with H being A, T or C) sites under different sequencing depths. The distribution of peaks in different genome components in each group, including the upstream 2 kb, 5'-untranslated region (UTR), 3'-UTR, CDS, intron, and downstream 2 kb was also analyzed.
The genomic differentially methylated regions (DMRs) between E50 and E60 groups were identified using the method described previously [49]. The number and the distribution of DMRs between two groups were analyzed. All DMR-containing genes were subsequently processed to functional enrichment analysis of GO and KEGG pathway using the DAVID (the Database for Annotation, Visualization and Integrated Discovery) web server. For GO and pathway analysis, Parent-child-Intersection method was used for enrichment analysis and Benjamini-Hochberg was used for multiple tests correction.

Bisulfite sequencing PCR (BSP)
There were three individual samples in E50 and E60 groups, respectively. Three relatively high methylated gene regions were chosen to validate MeDIP-seq data with BSP using individual samples. The primers for the three gene regions were listed in Additional file 1: Table S4. Bisulfite modification of 500 ng of genomic DNA was performed using Methyl Code Bisulfite Conversion Kit Fig. 6 Path-net and gene-gene interaction analysis of differentially methylated genes. a Path-net analysis of differentially methylated genes. Cycle nodes represent pathways, the size of the node represents the power of the interrelation among the pathways, and the lines represent the interaction between the pathways. The more edges in a pathway, the more pathways connecting to it and the more central role it plays within the network. b Gene-gene interaction analysis of differentially methylated genes. Red circle nodes represent up-regulated genes, and blue circle nodes represent down-regulated genes. The size of the node represents degree, and lines represent interactions between the genes. A straight line represents combinations anda dotted line represents indirect effects. a, activation; b, binding; c, compound; inh, inhibition; u, ubiquination; s, state change; dep, dephosphorylation