DNA methylation profiling of CpG islands in trigeminal ganglion of rats with orofacial pain induced by experimental tooth movement
Orofacial pain (OFP) refers to acute or chronic pain in orofacial regions such as teeth, maxillofacial soft tissue, alveolar bone, jaw, and temporomandibular joint. OFP has a high prevalence due to its multiple affected regions and thus brings serious medical consequences and social issues . Orthodontic treatment to correct malocclusion for better aesthetics and bite function is becoming increasingly well-known and popular. However, tooth movement-induced orofacial pain is the most cited negative effect during orthodontic treatment which hinders patients from seeking medical advice or adhering to treatment . Despite various treatment strategies such as nonsteroidal anti-inflammatory drugs (NSAIDs) have been applied to relieve orofacial pain due to tooth movement, none of them are truly effective without side effects . Given the high prevalence but limited treatment options of orofacial pain, research on the underlying mechanisms is crucial for developing more effective and secure treatment approaches. Numerous studies showed that the trigeminal ganglia (TG) plays a critical role in pain perception and its transmission in the orofacial regions. The trigeminal neurons situated in TG receive orofacial pain signals, which are then carried to the trigeminal nucleus caudalis. The signals are transmitted to the thalamus and are ultimately perceived by the sensory cortex of the brain . Additionally, a growing number of scientific literature demonstrated that DNA modification may play a crucial part in orofacial pain . DNA modification belongs to epigenetic regulation, mainly including DNA methylation and demethylation. In mammals, DNA methylation is a reversible and dynamic process regulated by DNA methyltransferases (DNMTs) and demethylases . Promoters, exons, and introns are the primary genetic locations that DNA methylation occurs. Particularly, CpG islands (CGIs) and CGI shores in promoter regions are the main occurrence sites, which often hinders transcription factor binding, thereby silences gene expression. It is reported that changes of DNA methylation status are involved in pathogenesis of several diseases . The differential expression of pain-related mediators (e.g., neuropeptides, cytokines, receptors, and ion channels) due to DNA methylation and demethylation in neurons, glial cells, and immune cells is instrumental in pain occurrence and development . Compared to studies on the role of DNA methylation in the dorsal root ganglion (DRG) and central nervous system (CNS), there are relatively few studies on DNA methylation related to in peripheral nervous system (PNS), especially the orofacial region and TG due to the more complex structure of the orofacial region. In complete Freund’s adjuvant (CFA)-induced orofacial inflammatory pain rats, the level of global DNA methylation, DNMT1 and DNMT3a expression in trigeminal ganglion were reduced due to peripheral inflammation. Additionally, some pain-related genes ( TRPA1 , TRPV1 , PIEZO2 , and P2X3 ) were confirmed to be epigenetic regulated through DNA methylation, while probably through DNMT3a in inflammation . Another research revealed significantly reduced DNA methylation of tumor necrosis factor-α ( TNF-α ) in the promoter region of TG in temporomandibular joint inflammatory pain, which resulted in the up-regulation of TNF-α . Despite the research on the role of some specific genes DNA methylation modulation, the DNA methylation profiling of CpG islands and CpG island shores in promoter regions in orofacial pain rats has not been elucidated yet. The aim of this study was to explore the methylation status of CGIs and CGI shores in promoter regions in TG of tooth movement induced orofacial pain rats, thus to further insight the DNA methylation regulation in orofacial pain. 

Animals and ethics Male Sprague Dawley rats (200–300 g) were purchased from Hangzhou Medical University. Each cage housed five rats, and rats were free access to water and food. The room for placing cages was temperature- and light-controlled with a 12-h:12-h light and dark cycle. The experimental procedures were approved by the Laboratory Animal Welfare and Ethics Committee of Zhejiang University (AP CODE: ZJU20220302). Establishment of tooth movement induced rat orofacial pain model A rat orofacial pain model was conducted by experimental tooth movement as reported previously . The rats first underwent general anesthesia via intraperitoneal injection of pentobarbital sodium (40 mg/kg). Subsequently, coil springs with a force of 40 g were ligated between the left incisor and the ipsilateral maxillary first molar. A dynamometer was applied to measure the force value. The rats in the Force group ligated by a coil spring with a force of 40 g, while the Sham group received the same coil spring without force. Pain assessment by rat grimace scale (RGS) Orofacial pain levels were evaluated by RGS according to previous studies . The rats were scored to obtain the baseline RGS scores before inserting the closed-coil springs. Briefly, the rats were placed separately in a transparent container (20.0 * 10.0 * 10.0 cm 3 ). After 15 min, the rats were observed and facial expressions were scored for 30 min (10 times, 3-min/time). This procedure was repeated 1, 3, 5, 7, and 14 days after closed-coil spring placement. The RGS assesses facial expression changes in four action units (eyes, ears, nose, and whiskers). Eye changes are characterized by orbital tightening; nose bulge is observed by a rounded extension of the skin on the nose bridge; ear position changes are ears pulled apart, angled forward or outward, or move back; whisker changes are characterized by backward, forward or clumped whisker position. A mean RGS score reflecting the pain level for each rat was obtained by averaging scores of the four action units of 10 measurements. The RGS score difference used in all analyses refers to the difference between the RGS scores after spring placement and baseline. Reduced representation bisulfite sequencing (RRBS) TGs of rats were isolated and collected from the 3 days Force and Sham groups. Genomic DNA was extracted following the manufacturer’s procedure by QIAamp Fast DNA Tissue Kit (Qiagen, Dusseldorf, Germany) and then sequenced by LC Technology (Hangzhou, China). The quality of the extracted genomic DNA was then measured and evaluated to ensure that only high-quality DNA was used for the following experiments. MspI (NEB, USA) was utilized to fragment DNA to approximately 200 bp samples, which were carried out to bisulfite conversion. Moreover, the Accel-NGS Methyl-Seq DNA Library Kit (Swift, MI, USA) was applied to connect adapters to single-stranded DNA fragments. In Brief, the steps of library construction and sequencing process mainly included the adaptase step, the extension step, and the indexing PCR step . Finally, sequencing was performed on Illumina Hiseq 4000 platform. Data analysis RRBS bioinformatic analysis was conducted by LC Technology (Hangzhou, China). Methylation statistics was performed to assess DNA methylation level, which was determined by the ratio of methylated C number to total number by in-house per scripts and MethPipe . R package-MethylKit was applied to calculate differentially methylated regions (DMRs) with a default parameter (1000 bp slide windows, 500 bp overlap, P < 0.05) . Moreover, the classification and function of differentially methylated CGIs and CGI shores annotated genes were detected using gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) analyses. To further focus on genes associated with pain, we matched and screened our differentially methylated CGIs and CGI shores annotated genes with the pain-related genes dataset obtained from GeneCards ( https://www.genecards.org ). Furthermore, the protein-protein interaction (PPI) network was achieved and established by Search Tool for the Retrieval of Interacting Genes (STRING) database to explore potential interactions of the selected annotated genes. The minimum required interaction score was set > 0.7 (high confidence) as the reliability value and the degree of connectivity of each protein was to assess its importance in the PPI network. Statistical analysis RGS data were expressed as mean ± standard error of the mean (SEM). Statistical analysis was performed using GraphPad Prism 8.0 software (GraphPad Software, San Diego, California USA). Unpaired t-test was used to compare between Force group and Sham group and P < 0.05 was considered statistically significant. 

Male Sprague Dawley rats (200–300 g) were purchased from Hangzhou Medical University. Each cage housed five rats, and rats were free access to water and food. The room for placing cages was temperature- and light-controlled with a 12-h:12-h light and dark cycle. The experimental procedures were approved by the Laboratory Animal Welfare and Ethics Committee of Zhejiang University (AP CODE: ZJU20220302). 

A rat orofacial pain model was conducted by experimental tooth movement as reported previously . The rats first underwent general anesthesia via intraperitoneal injection of pentobarbital sodium (40 mg/kg). Subsequently, coil springs with a force of 40 g were ligated between the left incisor and the ipsilateral maxillary first molar. A dynamometer was applied to measure the force value. The rats in the Force group ligated by a coil spring with a force of 40 g, while the Sham group received the same coil spring without force. 

Orofacial pain levels were evaluated by RGS according to previous studies . The rats were scored to obtain the baseline RGS scores before inserting the closed-coil springs. Briefly, the rats were placed separately in a transparent container (20.0 * 10.0 * 10.0 cm 3 ). After 15 min, the rats were observed and facial expressions were scored for 30 min (10 times, 3-min/time). This procedure was repeated 1, 3, 5, 7, and 14 days after closed-coil spring placement. The RGS assesses facial expression changes in four action units (eyes, ears, nose, and whiskers). Eye changes are characterized by orbital tightening; nose bulge is observed by a rounded extension of the skin on the nose bridge; ear position changes are ears pulled apart, angled forward or outward, or move back; whisker changes are characterized by backward, forward or clumped whisker position. A mean RGS score reflecting the pain level for each rat was obtained by averaging scores of the four action units of 10 measurements. The RGS score difference used in all analyses refers to the difference between the RGS scores after spring placement and baseline. 

TGs of rats were isolated and collected from the 3 days Force and Sham groups. Genomic DNA was extracted following the manufacturer’s procedure by QIAamp Fast DNA Tissue Kit (Qiagen, Dusseldorf, Germany) and then sequenced by LC Technology (Hangzhou, China). The quality of the extracted genomic DNA was then measured and evaluated to ensure that only high-quality DNA was used for the following experiments. MspI (NEB, USA) was utilized to fragment DNA to approximately 200 bp samples, which were carried out to bisulfite conversion. Moreover, the Accel-NGS Methyl-Seq DNA Library Kit (Swift, MI, USA) was applied to connect adapters to single-stranded DNA fragments. In Brief, the steps of library construction and sequencing process mainly included the adaptase step, the extension step, and the indexing PCR step . Finally, sequencing was performed on Illumina Hiseq 4000 platform. 

RRBS bioinformatic analysis was conducted by LC Technology (Hangzhou, China). Methylation statistics was performed to assess DNA methylation level, which was determined by the ratio of methylated C number to total number by in-house per scripts and MethPipe . R package-MethylKit was applied to calculate differentially methylated regions (DMRs) with a default parameter (1000 bp slide windows, 500 bp overlap, P < 0.05) . Moreover, the classification and function of differentially methylated CGIs and CGI shores annotated genes were detected using gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) analyses. To further focus on genes associated with pain, we matched and screened our differentially methylated CGIs and CGI shores annotated genes with the pain-related genes dataset obtained from GeneCards ( https://www.genecards.org ). Furthermore, the protein-protein interaction (PPI) network was achieved and established by Search Tool for the Retrieval of Interacting Genes (STRING) database to explore potential interactions of the selected annotated genes. The minimum required interaction score was set > 0.7 (high confidence) as the reliability value and the degree of connectivity of each protein was to assess its importance in the PPI network. 

RGS data were expressed as mean ± standard error of the mean (SEM). Statistical analysis was performed using GraphPad Prism 8.0 software (GraphPad Software, San Diego, California USA). Unpaired t-test was used to compare between Force group and Sham group and P < 0.05 was considered statistically significant. 

Tooth movement induced orofacial pain levels The rat orofacial pain model was induced by tooth movement with coil spring ligation between the left incisor and the ipsilateral first molar. Figure A demonstrated the rats at rest showing no orbital tightening, no nose bulge, ears standing straight, and whiskers in neutral position (black arrows). Figure B and C showed the rats in pain. The ears pull apart and back, whiskers cluster together and move backward (red arrows in Fig. B), as well as orbital tightening and obvious cheek bulge occur (red arrows in Fig. C). As displayed in Fig. D, the RGS score raised at 1 day, reached to peak at 3 days, then gradually decreased, and returned to the baseline at 14 days in both groups. However, the RGS score in the Force group was higher than that in the Sham group at all time points, with significant differences at 1 day (0.771 ± 0.085 vs. 0.404 ± 0.102, P = 0.02<0.05), 3 days (1.029 ± 0.124 vs. 0.408 ± 0.084, P = 0.002<0.01) and 5 days (0.563 ± 0.110 vs. 0.258 ± 0.044, P = 0.03<0.05). The most significant difference in orofacial pain was observed 3 days after ligating coil springs. Therefore, the isolated TGs were collected for RRBS after 3 days of coil springs ligation. Distribution of differentially methylated regions The distribution and number of total DMRs and DMRs in CGI and CGI shore in promoter region in tooth movement-induced orofacial pain rats compared with sham rats was analyzed and shown in Fig. . For total DMRs, 7022 promoter regions (19.43%), 16 exon regions (0.04%), and 29,104 intergenic regions (80.53%) were found hypermethylated regions, while 6857 promoter regions (18.29%), 11 exon regions (0.03%) and 30,616 intergenic regions (81.68%) were found hypomethylated regions (Fig. A). In addition, Fig. B showed that among 7022 hypermethylated promoter regions, 902 were promotor CGIs and 634 were CGI shores, while among 6857 hypomethylated promoter regions, 862 were CGIs and 577 were CGI shores. A volcano plot illustrated the filtered out DMRs in CGIs and CGI shores in promoter regions between the Force group and sham group, along with non-DMRs (Fig. A). As displayed in Fig. B, the heatmap further illustrated these DMRs. GO and KEGG enrichment analysis of differentially methylated CGIs and CGI shores in promoter regions GO analysis was performed to further study the DMRs of CGIs and CGI shores in the identified promoter regions. Figure A showed the GO analysis results, with the 3 top-ranked biological processes being negative regulation of transcription by RNA polymerase II, DNA strand renaturation, and angiogenesis. Neuron migration ranked 8th in the biological processes and might also be associated with pain. The 3 top-ranked molecular functions included nucleotide binding, RNA binding, and nucleic acid binding. The 3 top-ranked cellular components included the nucleus, cytoplasm, and membrane integral component. Figure B showed the top 20 GO terms in GO enrichment analysis, among which the 3 top-ranked terms were RNA polymerase II regulatory region sequence-specific DNA binding, negative regulation of transcription by RNA polymerase II, and retinoic acid binding. Furthermore, the function of differentially methylated CGIs and CGI shores in promoter regions was conducted using KEGG pathway enrichment analysis. Annotated gene symbols in the KEGG pathways that were found participated in several pain-related organismal systems including endocrine system, nervous system, immune system, and sensory system (Fig. C). Figure D displayed the top 20 pathways in enrichment analysis, among which 7 pathways (Fatty acid degradation, P = 0.00; Calcium signaling pathway, P = 0.01; Inositol phosphate metabolism, P = 0.02; Parathyroid hormone synthesis, secretion and action, P = 0.04; Wnt signaling pathway, P = 0.04; Neuroactive ligand-receptor interaction, P = 0.04; Cellular senescence, P = 0.04) were significantly enriched ( P < 0.05). Among the abovementioned 7 pathways, “Calcium signaling pathway”, “Wnt signaling pathway” and “Neuroactive ligand-receptor interaction” are usually considered to be involved in pain modulation. PPI analysis of pain-related genes After matching and screening our differentially methylated CGIs and CGI shores annotated genes with the pain-related genes dataset obtained from GeneCards ( https://www.genecards.org ), 730 genes (1651 methylated promoter CGIs and CGI shores) were identified to be related to pain mechanisms. A PPI network was established and revealed that several genes were important in the main clustering network, such as Ctnnb1 , Dlg4 , Creb1 , Camk2g , Bmp2 , etc. (Fig. ). Figure showed PPI network of associated pain-related genes enrichments in KEGG pathways including Calcium signaling pathway (Fig. A) and Wnt signaling pathway (Fig. B), and that in nervous system tissue expression (Fig. C). Among these genes, some showed hypo-methylated ( Ctnnb1 , Fgf2 , Hdac2 , Dlg4 , Creb1 , Bmp4 , Axin2 , Fzd4 , Wnt4 , Camk2g ) and others exhibited hyper-methylated ( Dvl1 , Creb1 , Vegfa , Myc , Fos , Igf1r , Wnt5a , Bmp2 , Bcl2l11 ), which may inhibit or promote the corresponding gene expression and contribute to tooth movement-induced orofacial pain. Table showed detailed information of different methylation status of pain-related example genes. According to the above findings, we speculated that DMRs in CGIs and CGI shores of promoter regions contribute significantly in the regulation of pain-related gene expression in TG of orofacial pain rats. 

The rat orofacial pain model was induced by tooth movement with coil spring ligation between the left incisor and the ipsilateral first molar. Figure A demonstrated the rats at rest showing no orbital tightening, no nose bulge, ears standing straight, and whiskers in neutral position (black arrows). Figure B and C showed the rats in pain. The ears pull apart and back, whiskers cluster together and move backward (red arrows in Fig. B), as well as orbital tightening and obvious cheek bulge occur (red arrows in Fig. C). As displayed in Fig. D, the RGS score raised at 1 day, reached to peak at 3 days, then gradually decreased, and returned to the baseline at 14 days in both groups. However, the RGS score in the Force group was higher than that in the Sham group at all time points, with significant differences at 1 day (0.771 ± 0.085 vs. 0.404 ± 0.102, P = 0.02<0.05), 3 days (1.029 ± 0.124 vs. 0.408 ± 0.084, P = 0.002<0.01) and 5 days (0.563 ± 0.110 vs. 0.258 ± 0.044, P = 0.03<0.05). The most significant difference in orofacial pain was observed 3 days after ligating coil springs. Therefore, the isolated TGs were collected for RRBS after 3 days of coil springs ligation. 

The distribution and number of total DMRs and DMRs in CGI and CGI shore in promoter region in tooth movement-induced orofacial pain rats compared with sham rats was analyzed and shown in Fig. . For total DMRs, 7022 promoter regions (19.43%), 16 exon regions (0.04%), and 29,104 intergenic regions (80.53%) were found hypermethylated regions, while 6857 promoter regions (18.29%), 11 exon regions (0.03%) and 30,616 intergenic regions (81.68%) were found hypomethylated regions (Fig. A). In addition, Fig. B showed that among 7022 hypermethylated promoter regions, 902 were promotor CGIs and 634 were CGI shores, while among 6857 hypomethylated promoter regions, 862 were CGIs and 577 were CGI shores. A volcano plot illustrated the filtered out DMRs in CGIs and CGI shores in promoter regions between the Force group and sham group, along with non-DMRs (Fig. A). As displayed in Fig. B, the heatmap further illustrated these DMRs. 

GO analysis was performed to further study the DMRs of CGIs and CGI shores in the identified promoter regions. Figure A showed the GO analysis results, with the 3 top-ranked biological processes being negative regulation of transcription by RNA polymerase II, DNA strand renaturation, and angiogenesis. Neuron migration ranked 8th in the biological processes and might also be associated with pain. The 3 top-ranked molecular functions included nucleotide binding, RNA binding, and nucleic acid binding. The 3 top-ranked cellular components included the nucleus, cytoplasm, and membrane integral component. Figure B showed the top 20 GO terms in GO enrichment analysis, among which the 3 top-ranked terms were RNA polymerase II regulatory region sequence-specific DNA binding, negative regulation of transcription by RNA polymerase II, and retinoic acid binding. Furthermore, the function of differentially methylated CGIs and CGI shores in promoter regions was conducted using KEGG pathway enrichment analysis. Annotated gene symbols in the KEGG pathways that were found participated in several pain-related organismal systems including endocrine system, nervous system, immune system, and sensory system (Fig. C). Figure D displayed the top 20 pathways in enrichment analysis, among which 7 pathways (Fatty acid degradation, P = 0.00; Calcium signaling pathway, P = 0.01; Inositol phosphate metabolism, P = 0.02; Parathyroid hormone synthesis, secretion and action, P = 0.04; Wnt signaling pathway, P = 0.04; Neuroactive ligand-receptor interaction, P = 0.04; Cellular senescence, P = 0.04) were significantly enriched ( P < 0.05). Among the abovementioned 7 pathways, “Calcium signaling pathway”, “Wnt signaling pathway” and “Neuroactive ligand-receptor interaction” are usually considered to be involved in pain modulation. 

After matching and screening our differentially methylated CGIs and CGI shores annotated genes with the pain-related genes dataset obtained from GeneCards ( https://www.genecards.org ), 730 genes (1651 methylated promoter CGIs and CGI shores) were identified to be related to pain mechanisms. A PPI network was established and revealed that several genes were important in the main clustering network, such as Ctnnb1 , Dlg4 , Creb1 , Camk2g , Bmp2 , etc. (Fig. ). Figure showed PPI network of associated pain-related genes enrichments in KEGG pathways including Calcium signaling pathway (Fig. A) and Wnt signaling pathway (Fig. B), and that in nervous system tissue expression (Fig. C). Among these genes, some showed hypo-methylated ( Ctnnb1 , Fgf2 , Hdac2 , Dlg4 , Creb1 , Bmp4 , Axin2 , Fzd4 , Wnt4 , Camk2g ) and others exhibited hyper-methylated ( Dvl1 , Creb1 , Vegfa , Myc , Fos , Igf1r , Wnt5a , Bmp2 , Bcl2l11 ), which may inhibit or promote the corresponding gene expression and contribute to tooth movement-induced orofacial pain. Table showed detailed information of different methylation status of pain-related example genes. According to the above findings, we speculated that DMRs in CGIs and CGI shores of promoter regions contribute significantly in the regulation of pain-related gene expression in TG of orofacial pain rats. 

As the earliest discovered epigenetic regulation modes, DNA methylation and demethylation are involved in pain regulation through altering expression of pain-related genes . Although abundant literature  indicates that altered DNA methylation contributes to pain sensitivity under various pain conditions, the majority of studies focused on DRG and the spinal cord rather than the orofacial region and TG. DNA methylation profiling of CpG islands and CGI shores in promoter regions in TG of the orofacial pain rat model has not been reported yet. Therefore, it was examined in the TG of experimental tooth movement-induced orofacial pain rats using RRBS. The RRBS results revealed that DMRs were mainly found and presented in intergenic regions and promoter regions, whereas exons and introns showed almost none, which may be attributed to the limitation of the RRBS technique. In comparison with whole-genome bisulfite sequencing (WGBS), RRBS is a technique that reduces the portion of the genome with high GC content by enzymatically digesting DNA, which allows for greater read depth. By restriction enzyme digestion, high-density CpG sites, such as CpG islands, can be sequenced with higher coverage. However, this approach examines a smaller fraction of the genome. Generally speaking, in RRBS, the total percentage of genome assessed < 20%, while the CpG/100 bp density > 3 . DNA methylation at different sites exerts various effects on gene expression. CGIs in the promoter region usually remain in a demethylated state, and gene transcription is usually hindered by methylation in the promoter region. In contrast, gene activity is increased in intragenic regions . However, DNA methylation does not only occur in CpG islands. The methylation and transcriptional activation of CpG island shore (i.e., the low CpG density region contiguous to the CpG island) are closely related . Therefore, this study focused on the DMRs in CGIs and CGI shores of promoter regions. In total, 2975 differentially methylated CGIs and CGI shores in promoter regions were identified. Among these, 902 CGIs and 634 CGI shores were found hypermethylated, and 862 CGIs and 577 CGI shores were found hypomethylated. GO and KEGG analyses were then performed to determine the distribution and function of annotated genes with differentially methylated CGIs and CGI shores in promoter regions, meanwhile, PPI network of pain-related genes was established to explore potential interactions of the selected annotated genes and discover the relatively important genes. As expected, GO analysis showed that the annotated genes participated in organismal systems including endocrine system, nervous system, immune system, and sensory system. In PPI analysis, nervous system was one of the gene enriched expression tissue, mainly including Creb1 , Dlg4 , Gad1 , Camk2g , Nos3 , Gnas , Fgf2 , Dnm1 , Fos , Gabbr1 , Gabrb3 , etc. KEGG pathway enrichment analysis revealed that “Calcium signaling pathway”, “ Wnt signaling pathway” and “Neuroactive ligand-receptor interaction” are significantly enriched pathways, which were previously reported to be related to various types of pain modulation . The pathways involved in orofacial pain were further verified by PPI clustering analysis. Specifically, functional enrichments in our network showed that Calcium signaling pathway mainly includes Camk2g , Adcy3 , Gnas , Cacna1s , Cacna1i , Adcy4 , Htr4 , Gna14 , etc. and Wnt signaling pathway mainly includes Ctnnb1 , Wnt5a , Wnt4 , Wnt9a , Axin2 , Dvl1 , Fzd1 , Fzd4 , Fzd9 , Myc , Lef1 , Ror2 , Sfrp2 , etc. Among the listed genes in Table , a special phenomenon was found that Creb1 exhibited both hypo-methylated and hyper-methylated on different gene position in promoter regions. Creb are transcription factors that stimulate the transcription of target genes whose docking has been shown to be inhibited by DNA methylation . Moreover, although opposing hypo-5-methylcytosine and hyper-5-hydroxymethylcytosine changes seem to co-localize mostly within gene bodies was reported previously , specific mechanisms remain to be explored. Among the aforementioned significantly enriched pathways and PPI network, several annotated genes were notable which might be associated with the causes or effects of tooth movement-induced pain due to the hypomethylation or hypermethylation of CGI or CGI shores in promoter regions of these genes. Ctnnb1 ( Catenin Beta 1 ) is key downstream part of the canonical Wnt signaling pathway. Ctnnb1 was proved to be one of the most enriched genes in neuropathic pain of brain by RNA-Seq . The Wnt/β-catenin pathway may be involved in pain by regulating cytokines. β-catenin was detected significantly raised in DRG of neuropathic pain model rats, while the pain was alleviated and TNF-α as well as IL-18 levels were decreased after using Wnt/β-catenin pathway inhibitor XAV939 . Wnt signaling pathways and some Wnt proteins are indispensable part of synaptic plasticity and pathogenesis of pain . Therefore, apart from Ctnnb1 , other genes such as Wnt5a , Wnt4 , Wnt9a , Axin2 , Dvl1 , Fzd1 , Fzd4 , Fzd9 , Myc , Lef1 , etc. may also play vital role in pain. Camk2g , Gnas , Cacna1s , Cacna1i , Adcy4 , Htr4 , etc. were genes enriched in Calcium signaling pathway. Camk2g , also known as calcium/calmodulin dependent protein kinase II gamma ( CaMKIIγ ), activates calmodulin-dependent protein kinase, shows activity in neuron projection, as well as participates in the regulation of nociception . Cacna1i ( Calcium Voltage-Gated Channel Subunit Alpha1 I ), encoding the Cav3.3 channel, is one of the main subtypes of T-type calcium channels that are widely expressed in the nervous system . T-type channels are considered potential targets for pain therapy, as abundant evidence showed that T-type channels participated in pathological pain sensation . Significantly increased Cacna1i expression was observed in TG isolated from foramen rotundum inflammatory constriction trigeminal infra orbital nerve mice . It was also reported as a causal gene of familial hemiplegic migraine . Gamma-aminobutyric acid A receptor beta 3 ( Gabrb3 ) is a pentameric ligand-gated chloride channel that mediates fast inhibitory signaling in neural circuits . Gabrb3 was reported to be necessary for the function integration of pyramidal neuron subtypes in the somatosensory cortex , and deletion of Gabrb3 in peripheral somatosensory neurons would cause mechanosensory dysfunction through loss of presynaptic inhibition of inputs to the central nervous system mediated by GABAA receptors . GABA and GABAA receptors play important roles in somatosensation and in the production of analgesic sensations . Htr4 is a serotonin receptor which can stimulate the production of cAMP in response to serotonin. Meanwhile, it modulates the release of various neurotransmitters in PNS and CNS . A recent study  showed decreased in Grm4 , Htr4 , Adrb2 , Kcnf1 , and Gad2 gene expressions and hypermethylation of these genes in chronic constriction injury (CCI) rats. Moreover, the study pointed out that 5‑Azacytidine may be a therapeutic target for alleviating pain sensation because it could increase expression of abovementioned genes through reducing their methylation levels. Our results preliminarily described the variation tendency of DNA methylation in tooth movement induced orofacial pain. However, limitations should also be noted. Firstly, the relatively small number of animals limits our study which may reduce the statistical power. Secondly, only young male rats were included which not allowed us to observe whether age and gender would influence the alteration of DNA methylation. Moreover, the current study is merely a preliminary report lacking functional verification of the observed candidate genes. Therefore, the detailed underlying mechanism of DNA methylation regulation involved in orofacial pain requires further investigation to subsequently validate and will be explored in future research, to yield more solid and innovative results, thus to provide an experimental basis for treatment of orofacial pain induced by tooth movement. 

Widespread methylation changes of CpG islands and CGI shores in promoter regions involving certain pathways and novel genes were found to be possibly related to orofacial pain induced by tooth movement, thus providing a basis for further study on the mechanism of DNA methylation in orofacial pain. 

