Whole-genome methylation profiling reveals regions associated with painful temporomandibular disorders and active recovery processes.

Data in this study were from participants in two protocols of the Orofacial Pain: Prospective Evaluation and Risk Assessment (OPPERA) study: a) a case-control study of chronic painful TMD (OPPERA I), and b) longitudinal study of people with recent onset painful TMD (OPPERA II). DNA specimens for the case-control study of chronic painful TMD were obtained from volunteers previously recruited for the OPPERA study,[ 40 ] and enrolled at four university research clinics in the eastern United States (Baltimore, MD; Buffalo, NY; Chapel Hill, NC; Gainesville FL) as previously described.[ 58 ; 60 ] Initiated in 2006, the OPPERA study is a population-based study of painful TMD including a prospective cohort of 3,263 participants (57% female, 43% male; 52% white, 32% African-American, 16% other race/ethnicity) aged 18–44 years who did not have a history of diagnosable TMD at enrollment. An additional 1,008 volunteers with examiner-verified painful TMD were simultaneously recruited at the same study sites, with pain symptoms experienced for at least 6 months. TMD was diagnosed using a modified version of the Research Diagnostic Criteria for TMD (RDC-TMD).[ 14 ] Information on other chronic pain conditions was collected by self-report on all participants using the Comprehensive Pain and Symptom Questionnaire (CPSQ).[ 45 ] We arrived at the number of samples needed to assay based on preliminary power calculations. From the full cohort, we first identified samples for methylation assay based on the available DNA concentration (>= 20 ng/μL) to only include samples that would be successfully measured. From the pool of samples with sufficient quantity, we randomly selected samples to minimize potential bias. All participants provided written, informed consent for all study procedures. The OPPERA parent study was reviewed and approved by institutional review boards at each of the study sites and at the data coordinating center, Battelle Memorial Institute. Further details of OPPERA study protocols[ 58 ] and primary findings[ 17 ] of the case-control study have been reported. Peripheral blood samples for longitudinal analysis were obtained from volunteers recruited for the OPPERA-II study.[ 55 ] These subjects were enrolled into a prospective, observational study of acute to chronic TMD. A telephone screening of 100,000 households in the counties surrounding the four OPPERA recruitment sites was used in a community-based sampling strategy. Subjects were enrolled as acute TMD cases if they reported orofacial pain for ≥ 5 days in at least one of the previous three months. Subjects were excluded if they experienced orofacial pain for more than 1 day/month in the 12 months preceding the current episode, or if they had any lifetime history of disabling orofacial pain. People identified as symptomatic during the telephone interview attended a three-hour clinic visit where TMD case-classification was determined by a trained examiner following Diagnostic Criteria for Temporomandibular Disorders (DC/TMD) criteria.[ 56 ] Participants meeting DC/TMD criteria were included in our analyses as recent onset TMD cases. All subjects also completed an extensive battery of psychosocial questionnaires and quantitative sensory testing, and each participant provided a whole blood sample at the same visit. Participants who completed the baseline clinic visit (V1) were asked to return for a follow-up visit six months after the baseline visit (V2), with a window of ±21 days. Assessment procedures were similar to baseline, including follow-up questionnaires completed prior to the visit, and most of the clinic procedures done at the baseline visit including the DC/TMD exam. A sample of whole blood was collected at the V2 as per baseline procedures. In our analysis, there is no participant overlap between OPPERA I samples and OPPERA II samples. TMD pain cases from OPPERA I cohort analyzed in this study all had chronic pain, i.e. < 6 months. This study was designed to detect patterns of methylation that distinguish chronic painful TMD cases from TMD-free controls. As an inclusive community-based sample intended to investigate risk factors including comorbid pain conditions, the OPPERA prospective study enrolled individuals both with and without concomitant pain. Therefore, we constructed two overlapping sets of reference groups in the analysis of chronic TMD, defining all TMD-free subjects as “non-TMD reference group” (non-TMD) and subjects free of TMD and any other chronic pain as “non-pain reference group” (non-pain). For the longitudinal study, recent onset cases (i.e. acute TMD cases at V1) were subsequently sub-classified according to V2 diagnosis: Persistent TMD cases if they had a positive TMD diagnosis at V2, or Transient TMD cases if they were negative for TMD at the follow-up visit. Demographic characteristics of the analytical groups were analyzed for group differences using Fisher’s exact test (for binary traits, e.g. sex), Student’s t-test (for nominal traits, e.g. age), or Chi-squared test (for categorical traits, e.g. race and smoking status). Peripheral whole blood was collected from all subjects by venipuncture into 5mL EDTA containing vacutainer tubes and stored at −80°C, as described previously.[ 59 ] Genomic DNA purification was performed using the Qiagen Extraction Kit (Qiagen, Hilden, Germany) following the manufacturer’s recommended protocol. DNA quality was measured using the Agilent Tapestation 4200 system (Agilent Technologies, Santa Clara, CA, USA). The Tapestation instrument provides a DNA integrity number, or DIN, 1 – 10 for each sample; any sample with a DIN of 7 or higher is considered to be of good quality and was subjected to bisulfite conversion. To accurately quantitate DNA we used Quant-iT Picogreen dye (ThermoFisher Scientific, Waltham, MA, USA) to stain the dsDNA, and then used the Tecan Infinite Pro-200 microplate reader (Tecan Group Ltd., Männedorf, Switzerland) to measure the florescence of the dyed dsDNA to produce an accurate quantitation. A total of 15μl of purified DNA were randomized across studies and case status and plated for whole genome methylation profiling using the Illumina Infinium MethylationEPIC platform (Illumina, San Diego, CA) by the Genome Quebec Innovation Center (Montreal, Canada). Bisulfite conversions were performed on ≥500ng of DNA per well using the EZ-96 DNA Methylation Kit (Zymo Research, Irvine, CA, USA), with the manufacturer’s recommended protocol optimized for downstream Infinium Methylation assays. For each sample, 4μl of bisulfite converted DNA was denatured and amplified overnight for 20 −24 hours. Fragmentation, precipitation and resuspension of the samples follow overnight incubation. After resuspension, samples were then hybridized to the Illumina Infinium Methylation EPIC BeadChip for 16–24 hours. Finally, the BeadChips were washed to remove any unhybridized DNA and then labeled with nucleotides to extend the primers to the DNA sample. Following the Infinium HD Methylation protocol, the BeadChips were imaged using the Illumina iScan system (Illumina, San Diego, CA, USA) and analyzed in Genomestudio v2011 with Methylation module version 1.9.0, and each CpG was mapped using GRCh37/hg19 genome assembly coordinates. Methylation signal data were checked for quality using Illumina’s Genome Studio software. Raw methylation signals, i.e. IDAT files, were imported with function read.metharray.exp from the Minfi package (version 1.36.0).[ 3 ] General quality control steps, including sex estimation, were conducted as implemented in the minfi package[ 3 ; 18 ] in R. Because our data was all collected from the same tissue (i.e., whole blood), we applied quantile normalization to methylation signals using minfi ’s preprocessQuantile function to minimize unwanted variation within and between samples.[ 66 ] As peripheral blood is a heterogeneous tissue, we estimated the relative proportion of pure cell types (NK-cells, CD4+ T-cells, CD8+ T-cells, monocytes, B-cells, and neutrophils) within each sample, following the method of Houseman[ 29 ] as implemented in minfi . The reference data set for peripheral blood composition was obtained from the FlowSorted.Blood.450k R package.[ 31 ] The matrix of cell composition estimates across all samples and cell types was retained as covariates for subsequent regression analyses. Methylation β-values were transformed (using a Logit function) to M-values for analysis; while β-values are more intuitive, they are subject to severe heteroscedasticity properties at highly methylated or unmethylated CpGs.[ 12 ] Linear regression was performed using the limma package in R (v. 3.46.0) to estimate the coefficient for the methylation difference between cases and the respective reference group at each CpG site (namely, DMC analysis), with adjustment for age, dummy-coded recruitment site, sex, smoking status, and cell composition as covariates. To account for racial stratification, we also included the first 3 principal components of the genotypes calculated previously using whole-genome SNP variant data[ 60 ] as covariates. Epigenome-wide significance was set at a conservative P < 5 × 10 −8 to correct for multiple testing.[ 53 ; 67 ] While differential methylation at individual CpG sites may influence gene expression, changes in methylation status across multiple consecutive CpG sites (differentially methylated regions, or DMRs) are more likely to reflect responses to environmental stimuli with biological consequences. We used the bumphunter[ 32 ] method in minfi to identify DMRs. This algorithm identifies regions of contiguous CpGs broadly associated with the trait of interest, which is more robust than identifying individual associated sites.[ 71 ] Candidate DMRs were examined for association by fitting a linear model adjusted for covariates as in the previous single CpG analysis. We used a bootstrap sampling method using 1000 bootstrap samples to compare the observed methylation differences at each DMR to a null model. Statistical significance was denoted by a family-wise error rate (FWER) < 0.1. Sex-stratified association tests were also performed as secondary analyses to investigate differences in methylation patterns between females and males. We performed pairwise DMR analyses in the recent onset TMD study as described above, to compare methylation across time points and between participants with persistent versus transient pain. For the primary analyses, we considered all DMRs detected by bumphunter and applied a conservative FWER < 0.1 to minimize false positives. As an exploratory analysis of significant DMRs from the analysis of chronic painful TMD, we used an unadjusted threshold of p < 0.05 to identify significant group differences in the comparisons between recent onset painful TMD cases and respective reference groups. To identify SNP variation likely to affect differential methylation of nearby CpGs (i.e., methylation quantitative trait loci, or mQTL[ 72 ]), we cross-referenced identified DMRs with a published dataset of replicated mQTLs. This analysis used data extracted from over 50,000 blood-based cis -mQTLs generated from the Brisbane Systems Genetics Study (BSGS) and the Lothian Birth Cohorts of 1921 and 1936.[ 42 ] We then extracted GWAS summary statistics for these SNPs (or proxies in high LD, r 2 ≥ 0.8) from two large cohorts characterized for orofacial pain to verify whether observed mQTLs were also potentially genetic risk factors for TMD. We used data from an analysis between facial pain cases (n=3,493), defined as participants who reported pain in the face or jaw for the past 3 months or more, and controls (n=163,843), defined as participants who did not report pain at any site in the past month.[ 39 ] The association test was constructed as a logistic regression model in REGENIE,[ 41 ] adjusted for age, age 2 , sex, first 40 PCs, and dummy-coded recruitment sites as covariates. UKB subjects were of white British ancestry, and included both males and females. We also looked for confirmation of mQTL association with TMD in a combined dataset incorporating four case-control studies of examiner-verified TMD: OPPERA, Sao Paulo Brazil TMD study (SPB), the OPPERA II chronic TMD replication study, and the Complex Persistent Pain Conditions (CPPC) study. [ 60 ] This analysis used a logistic regression model with a generalized mixed model to test for genome-wide association with TMD using the GMMAT package.[ 8 ] Age, sex, enrollment site and cohort were included as covariates, as well as the first 10 PCs, calculated on the merged raw genotypes, to account for population stratification. Differential methylation is a key epigenetic mechanism regulating gene expression, silencing or activating genes in response to environmental stimuli at critical expression quantitative trait methylation (eQTM) sites. We examined the effect of chronic painful TMD DMRs on the expression of nearby genes of a range of 200k basepairs downstream and upstream, in order to identify eQTM that may contribute to the development of TMD. Gene expression data was generated in the recent onset TMD (OPPERA II) cohort as described previously.[ 47 ] Genes were annotated to GRCh37/hg19 genome and featureCounts [ 36 ] was used to calculate gene expression. We applied transcript per million (TPM) values to measure and normalize gene expression levels. Quantile normalized M-values of CpG sites located in the chronic painful TMD DMRs were tested for correlation with the expression level of their overlapped and proximal genes. Association testing was performed with the lmFit function from the limma R package (version 3.46.0). A Benjamini and Hochberg adjustment was applied to correct for multiple testing and a false discovery rate (FDR) < 0.05 was adopted to screen for significant association. To examine whether methylation changes in blood are reflected in neuronal cells, we used the Blood Brain DNA Methylation Comparison Tool,[ 24 ] a publically available resource. This online database contains DNA methylation values for all probes present on the Illumina 450K Beadchip from 71–75 matched samples sourced from pre-mortem whole blood and four post-mortem brain regions (prefrontal cortex, entorhinal cortex, superior temporal gyrus, and cerebellum). One representative DMP from each of the six chronic painful TMD DMRs was extracted from the database.

A graphical overview of the study cohorts and analyses is presented in Fig. 1 . Demographic characteristics of these groups are reported in Table 1 . The final analytical dataset included 496 chronic TMD cases and 452 TMD-free controls, of which 342 were also considered pain-free controls as they reported no other pain conditions. The chronic TMD group included a higher proportion of females and non-White subjects than the non-TMD and non-pain groups. This group was also older and more likely to report a history of smoking than the reference groups. The acute TMD cohort included 62 recent onset TMD cases, of which 31 had Persistent TMD pain at V2 and 31 had no pain symptoms and were considered Transient TMD at V2. There were no significant differences in any of the demographic factors between the Persistent and Transient groups. Methylation analysis was performed on 1,140 specimens from the two OPPERA cohorts ( Fig. S1 ). Three samples were removed before preprocessing due to likely cross-contamination. Additional samples were removed for missing phenotypic data, mismatched sex, and belonging to excluded phenotype groups, leaving a total of 1107 samples for analysis ( Table S1 ). Finally, we removed 30,435 probes known to have common SNPs at the CpG site, 40,983 cross-reactive probes,[ 48 ] 17,920 probes from the sex chromosomes, and 25,495 probes that with low quality (detection p-value > 0.01) in at least one resulted in a total of 751,026 CpGs for analysis. We examined 751,026 individual CpG sites for association with chronic TMD case status. Separate analyses were conducted using different definitions for reference groups, either subjects without TMD (non-TMD cases) or subjects not reporting any pain conditions (non-pain cases) ( Fig. 2 ). A total of 187 differentially methylated CpGs (DMCs) were associated at a genome-wide threshold for significance of p < 5×10 −8 in the chronic painful TMD vs. non-TMD analysis. In the chronic painful TMD vs non-pain analysis, 701 DMCs exceeded the significance threshold. Of the 158 unique DMCs representing the intersection of these sets ( Table S2 ), 83 (52.5%) were hypermethylated in chronic painful TMD, while 75 (47.5%) had decreased methylation. For both analyses, we further noted enrichment of DMCs in compartments likely to have functional consequences for gene expression, such as within first exons or near transcription start sites, and depletion of DMCs in “open sea” regions ( Table S3 ). To identify functional patterns of methylation differences more likely to have robust effects on gene expression,[ 43 ] we implemented the bumphunter algorithm to find DMRs. The algorithm evaluated 600 regions as potential DMRs (217 regions from chronic painful TMD vs. non-TMD, 383 from chronic painful TMD vs. non-pain), of which six had an FWER < 0.10 in either the non-TMD or non-pain analysis ( Table 2 ). Two loci on chromosome 1 were associated with chronic painful TMD: one near the FMOD gene (Chr1:203320190–203320732; p=1.71×10 −3 ; Fig. 3a ), and another downstream of the PM20D1 gene (Chr1:205818956–205819609; p=2.12×10 −4 ; Fig. 3b ). One associated DMR was located on chromosome 4 (Chr4:124232–125504; p=6.09×10 −5 ; Fig. 3c ) near the ZNF718 gene locus. Two loci were located on chromosome 6, one in the region of the ZFP57 and MOG genes (Chr6:29648225–29649024; p=2.44×10 −3 ; Fig. 3d ), and the other near the RNF39 locus (Chr6:30039132–30039801; p=1.61×10 −3 ; Fig. 3e ). The sixth DMR was located in a region on chromosome 9 near the TRAF2 gene locus (Chr9:139781209–139781630; p=1.92×10 −3 ; Fig. 3f ). Separate sex-stratified analyses distinguished all of these regions as being associated with chronic painful TMD in either females only or males only, except for the DMR on chromosome 9 near TRAF2 ( Table S4 ). One additional locus, on chromosome 11 (Chr11:5617367– 5618408, p=8.91×10 −4 ) near the TRIM34 gene, was only found significant from analysis in females between chronic painful TMD and non-TMD, indicating the region was potentially specific to female chronic TMD patients. To understand the temporal dynamics of methylation changes in acute painful TMD, we performed a DMR analysis in the OPPERA II cohort of recent onset TMD, investigating pairwise contrasts between Persistent and Transient pain groups at V1 and V2, and within the Persistent and Transient groups between V1 and V2 time points. Ten DMRs were significant after correcting for FWER < 0.10 across these contrasts, including eight distinct loci ( Table 3 ). The majority of differences were observed at V2 between those with persistent pain versus those with transient pain, suggesting that most dynamic changes in methylation distinguishing individuals that improve their pain from those that do not occur during the course of pain becoming chronic. We then examined the six DMRs associated with chronic painful TMD for dynamic changes in methylation in this cohort, both to replicate patterns of methylation observed in chronic TMD and also to explore changes in methylation of these genes across the transition from healthy to chronic disease states ( Table 4 ). Overall, we found four of the six DMRs exhibited significant ( P < 0.05) methylation changes in recent onset TMD groups. We further inspected the direction of methylation changes of the four DMRs in the analyses of OPPERA I samples (comparisons between chronic painful TMD and non-TMD/non-pain) and OPPERA II samples (comparisons between persistent painful TMD and transient painful TMD at both visits). The four DMRs were differentially methylated in the same direction in the Persistent group at V1 as in OPPERA I analyses. At V2, only two of the four DMRs showed the same direction in the Persistent pain group. As lower levels of methylation are generally associated with the activation of gene transcription, these results predict higher levels of expression of these genes in individuals at risk of chronicity. This pattern corresponds to the direction of methylation in chronic painful TMD for three regions (the FMOD , ZNF718 , and ZFP57/MOG loci), in which TMD cases were associated with hypomethylation. However, the direction of differential methylation was the opposite in the PM20D1 region, in which hypermethylation was associated with chronic painful TMD caseness. Similarly, we examined the eight DMRs associated with acute painful TMD for their methylation status in the chronic painful TMD cohort compared to the TMD-free and pain-free reference groups ( Table 3 ). Two DMRs ( POU5F1 and CYP2E1 ) exhibited differential methylation (P < 0.05) in chronic painful TMD in comparisons with non-TMD or non-pain and this differential methylation was in the same direction as in recent onset painful TMD; they were hypomethylated in both chronic TMD and during the acute phase of TMD in those who continued to have persistent pain. These results suggest hypomethylation changes in the 6 months following an acute TMD episode, corresponding to the higher expression of respective genes, represents a risk factor for chronicity that is maintained in the chronic state. Genetic variation can influence CpG methylation, thereby exerting an indirect effect on gene expression beyond direct effects on promoter or transcription factor binding sites. We examined if CpG methylation in the chronic painful TMD DMRs was affected by nearby SNPs, as potential cis -acting mQTLs. We tested associated SNP-CpG pairs in the McRae et al. dataset,[ 42 ] querying all CpGs located within the 6 chronic TMD DMRs. A total of 53 SNP-CpG pairs were extracted ( Table S5 ), which represented 13 distinct mQTLs. These included 2 mQTLs for the FMOD region (rs7543148, rs10920616), 1 for the PM20D1 region (rs823080), 1 for the ZNF718 region (rs147169138), 3 for the ZFP57 region (rs374317, rs416568, rs2747421), and 6 for the RNF39 region (rs114942112, rs115407701, rs114141756, rs115226254, rs115119809, rs116518960). We next examined these 13 mQTLs, or nearby proxy SNPs in strong LD, for association with a TMD pain-related phenotype in two large cohorts, the UK Biobank and a mega-analysis of four examiner-verified TMD cohorts ( Table 5 ) to confirm whether mQTLs had observable functional consequence. In the UK Biobank, a proxy SNP representing 3 RNF39 region mQTLs was associated with orofacial pain (rs115226254, rs115407701, rs116518960, OR = 0.91, p = 0.04). In the mega-analysis of chronic TMD, a proxy SNP representing a PM20D1 region was associated with TMD (rs823080, OR = 1.14, p = 0.0041). DNA methylation is a mechanism for dynamically regulating gene expression that is sensitive to environmental stimuli and homeostatic pressures. The relationship between methylation changes and activation or repression of transcription within a cell type may be dependent on moderating factors specific to temporal or etiological states. We sought to understand how differential methylation dynamically regulates gene expression under different conditions. As the OPPERA II cohort provided specimens at two time points (i.e., recent TMD onset and 6 months later) and included groups with different outcomes (i.e. Persistent and Transient), we examined associations between methylation of CpG sites within our significant DMRs and mRNA expression of genes within 200 kb to identify dynamic eQTMs in cis . Expression-methylation relationships were examined in all available subjects during the acute painful TMD period (V1) and separately at 6 months later (V2). We also identified eQTMs longitudinally, combining both time points within Persistent and Transient TMD groups, respectively. This analysis confirmed that differential methylation of a particular DNA region has functional consequences for expression of nearby genes, often affecting more than just one neighboring gene, but that the effect of methylation is not always consistent between groups ( Table 6 ). Co-regulation between the respective DMR CpGs and the nearest gene was strong in most or all subgroups for cg07533224 (with PM20D1 ) and cg15570656 (with ZFP57 ) indicating a tight inverse relationship between methylation and expression of the expected, most proximal gene. Other CpGs were only associated with the nearest gene in particular subsets of the cohort. A correlation between cg26987645 and FMOD was only observed in Transients, possibly suggesting the relationship between methylation and FMOD expression is only uncovered in acute TMD cases engaging in active recovery processes. Similarly, ZNF595/ ZNF718 also only exhibits this correlation in the Transients group. Several CpGs are correlated with the transcript levels of genes in addition to, or more strongly than, the nearest gene, including cg26987645 with BTG2 and OPTC ; cg15570656 with MOG ; and cg13918754 with ZNRD1 . This may suggest regulation of these genes is a more critical function of the respective DMR. The CpG on chr9, cg25078751, was not strongly associated with expression of any nearby genes. Absolute levels of methylation extracted from the Blood Brain DNA Methylation Comparison Tool were similar between blood and brain regions for representative CpG sites in four out of the six DMRs, with moderate to high levels of correlation in five of the regions ( Fig. S2 ), suggesting that blood may be an expedient proxy tissue for tracking changes in neurons and other cells involved in pain processing. We note, however, that the Tool uses postmortem brain tissue, which may impact methylation profiles. For the FMOD region, methylation was much higher in the brain tissue than in blood, and the correlation was only moderate (r=0.48) but significant (P=1.77e-05). For the PM20D1 region, methylation of CpGs was highly variable between individuals in both blood and brain, but correlations between tissues was very strong (r=0.98). Methylation in the ZNF718 region tended to be very low for brain, and was similarly low in blood but more variability was observed; correlation between blood and brain was strong (r=0.811). In the ZFP57 / MOG region, most individuals had high levels of methylation in both blood and brain, but significant variability was observed with high correlation between tissues (r=0.93). Differences were observed in the RNF39 region, with lower mean methylation in blood relative to brain tissues; correlation between blood and brain was moderately high (r=0.77). For the TRAF2 region, methylation levels were low for all tissues examined, with no correlation between brain and blood (r=0.02).

Dynamic regulation of genes in response to environmental exposures (e.g., injury, trauma, psychosocial stress) is a critical mechanism underlying the long-term pathological changes that lead to chronic pain.[ 33 ] The primary findings of this study, the first to assess genome-wide DNA methylation differences associated with chronic TMD, include a number of genomic regions that may contribute to TMD etiology. Further analyses in a separate cohort of recent onset TMD revealed patterns of differential methylation that may determine persistence or resolution of orofacial pain. As regions of consistently differential methylation are more likely to be functional than single CpG sites,[ 25 ; 30 ] we focused our follow-up analyses on the 6 DMRs associated with chronic painful TMD and their neighboring genes, integrating the methylation datasets with other -omic data exploring genetic variation and transcriptomic features. These analyses provided strong convergent evidence to support the physiological significance and regulatory function of these DMRs. Several of the DMRs, including ZFP57 , RNF39 , and ZNF718 , likely modulate a variety of genes downstream as regulators of transcription or methylation. For example, individuals with chronic TMD exhibited hypomethylation at the promoter and first exon of the zinc finger protein 57 ( ZFP57 ) gene in the Major Histocompatibility Complex (MHC) region, as did individuals with persistent TMD pain relative to those with transient pain at both time points. This gene functions as a transcriptional repressor by maintaining methylation of other genes, especially during development as part of genomic imprinting. It may precipitate a cascade of methylation events in response to environmental stress or inflammation, as it has been described as an epigenetic mark of post-traumatic stress disorder (PTSD)[ 52 ; 68 ] and periodontitis.[ 26 ] A second DMR in the MHC region was located in the last intron and terminal exon of the RING finger protein 39 ( RNF39 ) gene. Hypomethylation in this gene locus was associated with chronic but not acute painful TMD at either visit. RNF39 is an immediate early gene involved in plasticity in the hippocampus, and participates in immune and inflammatory functions. Differential methylation at this locus has been associated with prenatal chemical exposure[ 61 ] and early-life experience of famine,[ 35 ] as well as with PTSD.[ 52 ] Several reports to date describe an association between the RNF39 locus and chronic pain, either hypomethylation with musculoskeletal pain,[ 44 ] or conversely, hypermethylation with chronic widespread pain following motor vehicle collision[ 6 ] and migraine.[ 21 ] We also observed hypomethylation in chronic painful TMD as well as in those individuals with acute TMD with persistent pain at 6-months post-onset at V1 for a DMR spanning the promoter and first exon of the zinc finger 718 ( ZNF718 ) gene. However, at V2, the direction for methylation of this region was the opposite. This gene is most likely a transcription factor and may function as a modulator of inflammatory response as it has been observed to be hypermethylated in periodontitis.[ 26 ] Biphasic regulation of this gene may be a reflection of its critical function in inflammatory response as we have recently described a biphasic pattern in the contribution of inflammation to pain resolution.[ 47 ] Thus, early hypomethylation of ZNF718 and consequent increase of its expression would reduce the risk for chronic pain; however, subsequent hypomethylation of the gene may be a marker of a pain state. The DMR located in the promoter region of the PM20D1 , or peptidase M20 domain containing 1, gene likely has a more complicated relationship with the inflammatory response. Hypermethylation was observed in chronic painful TMD and subjects at V1 with acute painful TMD with pain persisting at 6 months after onset. This was in agreement with studies reporting increased methylation associated with chronic musculoskeletal pain in adults[ 44 ] and post-surgical pain in children.[ 10 ] However, these findings are contrasted by hypomethylation at V2 in those with pain persisting at 6 months. Again, we observe biphasic regulation of an inflammatory gene. Establishing true causal relationships with epigenetic markers remains challenging, as cross-sectional studies of potentially dynamic methylation changes are unable to determine whether differences between cases and reference groups signify a cause or consequence of the pathophysiology underlying TMD. When we compared the methylation status of sequentially sampled recent-onset TMD participants at the methylome-wide level, we found that differential methylation was most strongly observed when comparing those whose pain persisted versus those whose pain was transient at the second visit V2 ( Table 3 ), supporting the hypothesis that changes to methylation status in response to acute pathology are largely associated with pain resolution. Two out of eight DMRs observed in the longitudinal study also exhibited methylation changes in chronic painful TMD, in the same direction. We thus were able to replicate two of our DMRs, at the CYP2E1 and POU5F1 loci, between two TMD cohorts and established the timeline of this differential methylation. These two genes, cytochrome P450 family member CYP2E1 responsible for drug metabolism and lipid synthesis, and POU5F1 , a transcription factor important in embryonic development and stem cell pluripotency, would require further investigation regarding their role in TMD pathogenesis. It is possible that the rest of the longitudinal DMRs represent transiently methylated regions critical for association with chronic painful TMD but not its maintenance. Conversely, four of the six chronic painful TMD DMRs also exhibited a significant change in methylation status in individuals with persistent pain relative to transient pain after recent onset ( Table 4 ); this would indicate stably regulated epigenomic sites contributing to activation of gene expression leading to chronic pain. These sites were already hypomethylated (for ZFP57 , RNF39 , and ZNF718 loci) or hypermethylated (for PM20D1 ) in the individuals with recent onset pain (V1 visit), and thus can serve as predictors of the development of chronic painful TMD. It remains unknown whether these sites were hypomethylated prior to the acute TMD episode and thus predispose for chronic pain, or if they were rapidly hypomethylated following TMD onset. We further characterized the functionality of the identified chronic painful TMD DMRs to understand their impact on etiological processes. We first identified 13 mQTLs, representing genetic variations at the SNP level associated with differential methylation at five of the six chronic painful TMD DMRs. Furthermore, we found corroborating evidence for two of these mQTLs in other cohorts for association with orofacial pain. As the vast majority of disease-associated SNPs are not located in coding regions of genes,[ 57 ] it is likely that many of them exert their effects by modulating gene expression via methylation. Our findings support this supposition, as we also established a strong inverse relationship between the majority of the chronic painful TMD DMRs (all except the DMR on Chr 9) and gene expression of one or more nearby genes. Our results should inform the design of future studies of differential methylation and pain in a heterogeneous convenience tissue. Blood was used in this study as a convenient and informative proxy source of DNA, but our findings should largely be representative of methylation changes in other tissues, as supported by the generally strong correlations between blood and brain tissue methylation observed for our DMRs. As such, these DMRs may be useful as measurable biomarkers of disease processes. Regardless, the best tissue type to assess for TMD pathology is not known, and given the complexity of the phenotype, there may be no single best tissue. The power of our chronic painful TMD case-control analysis was sufficient to implicate DMCs and DMRs at genome-wide significance for the expected effect size; however, the observed differences in methylation specific to TMD in blood samples were typically smaller than often assumed, so sample sizes for chronic pain EWAS may need to be higher than other traits or conditions. We observed differences in β-values of 2–3% in our significant DMRs, whereas power calculations for EWAS typically assume higher changes (8–10%). The size (n=62) of the recent onset of painful TMD cohort also limited our power in the gene expression analysis and other potential integrative analyses. We also chose to use reference datasets for detection of mQTLs due to potentially lower power for discovery in our dataset alone, but using such external sources may affect the results as the analytical parameters (e.g. covariates) of these reference datasets differ from our data. Finally, the generalizability of our findings is limited by the restricted age range and geographic/demographic characteristics of our study cohorts. Likewise, we recognize that there are assumptions and limitations built into the methylation assay and bumphunter algorithm used to identify DMRs that may affect our results. However, our use of a widely utilized platform and the deposition of our data into a public database will facilitate potential replication of our results by other groups. In summary, this study demonstrated that differential DNA methylation at a number of genetic loci is associated with painful TMD, with evidence for dynamic regulation over the course of onset and development of the disorder. Our preliminary functional studies support an etiological role for these DMRs in activating inflammatory pathways that initiate and alternately maintain or resolve chronic pain. These findings support further investigation into strategies to interrupt these etiological pathways at critical points in the development of chronic pain following trauma or injury.

Temporomandibular disorders (TMD) are the leading type of non-odontogenic orofacial pain and are frequently comorbid with other chronic pain conditions, including fibromyalgia[ 1 ; 34 ], headache, back pain, and stomach pain, [ 13 ; 27 ] suggesting common pathophysiology. Like other types of chronic pain, it is believed that both genetic and environmental factors affect the susceptibility to TMD; the heritability of TMD has been estimated between 17–27%[ 49 ; 60 ] in family and community-based studies, with the discovery of a number of associated genes.[ 54 ; 60 ] Despite the increasing evidence for genetic factors, to date the observed genetic risk markers explain only a tiny portion of genetic variance estimated in the population. Epigenetic factors may be an alternate mechanism contributing to variability in TMD susceptibility. Tissue-specific expression of genes is dynamically regulated by methylation of cytosine at CpG dinucleotides.[ 20 ] DNA methylation controls the binding of regulatory proteins such as transcription factors and silencing complexes. DNA methylation in gene promoter regions is generally associated with a decrease in translation, although the functional consequences of such modification are complex. Across tissue types, methylation is dynamically regulated in response to environmental exposures such as injury, toxins, and infection, suggesting that methylation may be a critical component of pro-nociceptive states.[ 19 ] With the recognition that DNA methylation is variable in adult, non-dividing cells,[ 70 ] such changes in methylation are accompanied by differences in pain behavior.[ 19 ; 37 ; 50 ; 64 ; 65 ] In animal models, methylation differences in a number of genes have been associated with specific pain states, including cystathione-B-synthetase (CBS),[ 51 ] miR-219,[ 46 ] and TNFα[ 4 ] in inflammatory pain; P2rx3,[ 73 ] Wnt3a,[ 16 ] and Gad1[ 69 ] in neuropathic pain; and Nr3c1, Cnr1, and Trpv1 in visceral pain.[ 28 ] With limited access to neuronal and other local tissues in human pain patients, candidate gene and epigenome-wide association studies (EWAS) using whole blood samples have identified methylation differences associated with chronic pain conditions including endometriosis,[ 15 ] fibromyalgia,[ 11 ] cancer pain,[ 63 ] widespread musculoskeletal pain,[ 38 ] migraine,[ 21 ] low back pain,[ 2 ; 23 ] and complex regional pain syndrome.[ 7 ] Longitudinal changes in DNA methylation are beginning to be assessed in patients undergoing surgery, although larger patient samples are needed.[ 9 ; 10 ] Differential methylation of the promotor region of the TRPA1 gene, which encodes a temperature-sensing ion channel, is associated with both heat pain sensitivity[ 5 ] and pressure pain threshold,[ 22 ] suggesting that methylation-based gene expression may contribute to pro-nociceptive states leading to the development of chronic pain. Methylation patterns may serve as biomarkers of distinct etiological pathways, distinguishing patients with different underlying pain mechanisms.[ 62 ] Despite these findings, advances in understanding the role of DNA methylation in pain has been hindered by the lack of sufficiently powered studies, and by the imperfect correlation between methylation states in neurons and those observable in available blood samples.[ 24 ] In the present study, we performed an EWAS to detect sites and regions of differential methylation between individuals with and without TMD using whole peripheral blood, in order to identify assessable epigenetic biomarkers of both acute and chronic orofacial pain. We then explored the functional consequences of genomic variation and methylation-based regulation on gene expression of implicated genes.