Genome-wide association reveals contribution of MRAS to painful temporomandibular disorder in males.

Painful temporomandibular disorders (TMDs) are the leading cause of chronic orofacial pain, but its underlying molecular mechanisms remain obscure. Although many environmental factors have been associated with higher risk of developing painful TMD, family and twin studies support a heritable genetic component as well. We performed a genome-wide association study assuming an additive genetic model of TMD in a discovery cohort of 999 cases and 2031 TMD-free controls from the Orofacial Pain: Prospective Evaluation and Risk Assessment (OPPERA) study. Using logistic models adjusted for sex, age, enrollment site, and race, we identified 3 distinct loci that were significant in combined or sex-segregated analyses. A single-nucleotide polymorphism on chromosome 3 (rs13078961) was significantly associated with TMD in males only (odds ratio = 2.9, 95% confidence interval: 2.02-4.27, P = 2.2 × 10). This association was nominally replicated in a meta-analysis of 7 independent orofacial pain cohorts including 160,194 participants (odds ratio = 1.16, 95% confidence interval: 1.0-1.35, P = 2.3 × 10). Functional analysis in human dorsal root ganglia and blood indicated this variant is an expression quantitative trait locus, with the minor allele associated with decreased expression of the nearby muscle RAS oncogene homolog (MRAS) gene (beta = -0.51, P = 2.43 × 10). Male mice, but not female mice, with a null mutation of Mras displayed persistent mechanical allodynia in a model of inflammatory pain. Genetic and behavioral evidence support a novel mechanism by which genetically determined MRAS expression moderates the resiliency to chronic pain. This effect is male-specific and may contribute to the lower rates of painful TMD in men.

1. Introduction

Painful temporomandibular disorder (TMD) is the most common cause of chronic orofacial pain[8] and is characterized by pain (including allodynia and hyperalgesia), jaw dysfunction, fatigue, and psychological disturbance.[7] Temporomandibular disorders affect about 5% of adults in the United States,[13,29] with an annual incidence of 3.9% for pain-related TMD.[38] A broad range of factors predict higher risk of developing painful TMD, including: greater age, presence of comorbid pain conditions,[1] greater frequency of somatic symptoms,[10] and poor sleep quality.[34]

Evidence supports a role for genetic contributions to TMD, with twin studies demonstrating heritability of 27%.[30] In a case–control study using a candidate panel of 350 genes relevant to nociceptive, inflammatory, immunological, and affective processes, we found that TMD was associated with single-nucleotide polymorphisms (SNPs) in 7 genes: COMT, an enzyme catabolizing catecholamines; HTR2A, one of the receptors for serotonin; NR3C1, the glucocorticoid receptor; CAMK4, the calcium/calmodulin-dependent protein kinase; CHRM2, the muscarinic acetylcholine receptor; IFRD1, the interferon-related developmental regulator; and GRK5, the G-protein-coupled receptor kinase.[41] To date, only one genome-wide association study (GWAS) of orofacial pain has been published.[35] Using self-reported orofacial pain as the phenotype, it found 2 female-specific genome-wide significant associations that were subsequently replicated in independent cohorts. One SNP was near the gene encoding relaxin/insulin-like family peptide receptor 2, RXP2; the second SNP was in the BAH domain and coiled-coil–containing gene, BAHCC1.

The aims of this study were to undertake the first discovery GWAS of examiner-verified chronic TMD, evaluate generalizability of the findings through independent replication, and investigate the functional biology of genetic markers that emerged from discovery and replication phases. The discovery sample included 3104 racially diverse participants from the Orofacial Pain: Prospective Evaluation and Risk Assessment (OPPERA) study, where painful TMD was determined by clinical examination using modified Research Diagnostic Criteria for TMD (RDC/TMD), and genotyping was conducted using the Illumina HumanOmni 2.5M bead chip. Single-nucleotide polymorphisms that achieved genome-wide significance (P < 5 × 10−8) in the discovery cohort were carried forward for replication in 7 additional cohorts. Replication studies included TMD case–control studies as well as population-based, cross-sectional studies, and were combined by meta-analysis to report size and direction of genetic effects. We then investigated the functional consequences of associated SNPs on etiological pathways using bioinformatic and transgenic approaches.

2. Methods

An overview of the study design, including discovery, replication, and validation phases, is provided in Figure 1.

Figure 1.

Study design overview diagram. GWAS, genome-wide association study; KO, knockout.

2.1. Discovery genome-wide association study cohort

The discovery case–control cohort included 3104 subjects selected for genome-wide genotyping, drawn from the OPPERA study (n = 2868) and a separate case–control study of chronic TMD (n = 236). The OPPERA subjects were men and women (61.9% women) between the ages of 18 to 44 years, recruited at 4 sites in the eastern United States (Buffalo, NY; Baltimore, MD; Chapel Hill, NC; and Gainesville, FL); enrollment of this cohort has been described previously.[37] The OPPERA subjects were racially and ethnically diverse: 57.5% white, 25.7% black, and 16.8% other/not reported (Supplementary Table 1, available at http://links.lww.com/PAIN/A688). From the full OPPERA cohort, subjects selected for GWAS analysis were determined randomly, leaving approximately 1000 subjects for later replication studies. The additional 236 subjects were from a separate case–control study of chronic TMD conducted by the same OPPERA study team.[39] Subjects were all non-Hispanic white women, aged 18 to 44 years, recruited at the Chapel Hill, NC (UNC) site. The combined group reported in this article consisted of 1011 TMD cases and 2093 TMD-free controls. All participants provided informed, signed consent for all study procedures. The OPPERA study was reviewed and approved by institutional review boards at each of the 4 study sites and at the data coordinating center, Battelle Memorial Institute. The TMD case–control study was approved by the Biomedical Institutional Review Board of the University of North Carolina at Chapel Hill.

2.2. Case classification

The classification of TMD[27] was based on the RDC/TMD.[7] Cases met all 3 of the following criteria: (1) pain in the cheeks, jaw muscles, temples, or jaw joints that occurred for at least 5 days per month during the preceding 6 months, including at least 15 days in the month before enrollment; (2) pain reported in the examiner-defined orofacial region for at least 5 days out of the prior 30 days; and (c) pain evoked by palpation of the orofacial muscles or maneuver of the jaw that occurred in at least 3 masticatory muscles or at least 1 temporomandibular joint or both. Examiners were trained in standardized procedures used during the examination, and annual evaluation of interexaminer reliability revealed consistently excellent levels of examiner agreement (kappa > 0.75).

2.3. Genotyping and imputation

Study participants provided a sample of whole blood that was stored in 5 mL EDTA-containing polyethylene vacutainers. Genomic DNA purification was performed using the Qiagen Extraction Kit following the manufacturer's recommended protocol. Genotyping was performed by the Center for Inherited Disease Research using the Illumina HumanOmni 2.5 Exome Bead Chip platform (Illumina, Inc, San Diego, CA). Genetic data cleaning was accomplished by the Genetic Analysis Center at the University of Washington following their established pipeline.[19] Batch effects were assessed by comparing missing call rate and the χ2 test for allelic frequency between genotyping batches. Sample identity and sample quality analyses included checks for missing call rate, chromosomal anomalies, cryptic relatedness, autosomal heterozygosity outliers, gender mismatch, and genetic ancestry. Consistency of genotyping calls was assessed using 69 duplicates of study samples and 66 stock samples from HapMap reference subjects. Single-nucleotide polymorphism level quality checks included assessments of missing call rate, duplicate discordance, and Mendelian errors. Because of the mixed population structure of the OPPERA cohort, Hardy–Weinberg equilibrium was tested in homogenous European- and African-ancestry groups separately.

Cleaned genotypes were then imputed to the 1000 Genomes Project phase 3 reference panel using the software packages SHAPEIT[6] for prephasing and IMPUTE2[11] for imputation. Instead of using the quality metric provided by IMPUTE2, imputation quality was assessed by calculating the OEvar statistic, which is the ratio of the observed variance of imputed dosages to the expected binomial variance,[21] and which ranges from 0 to 1. Imputed SNPs were excluded from the analysis based on a filter for effectively small sample size, taking into account the number of minor alleles as well as imputation quality. This filter was calculated using the formula: 2 × MAF × (1 − MAF) × N × OEvar ≥ 50, where N is the number of TMD cases in the cohort (ie, 999), MAF is the minor allele frequency, and the threshold for inclusion was set at 50 for binary traits.

2.4. Data availability

Study data have been deposited and made publicly available at the Database of Genotypes and Phenotypes (dbGaP) public repository (accession number phs000796.v1.p1).

2.5. Genotype quality assessment

The discovery data set included the maximal set of unrelated individuals passing all quality checks, and consisted of 999 TMD cases and 2031 TMD-free controls. Initial genotyping results were provided for 2.57 million SNPs; a total of 2.53 million SNPs passed quality filters including call rate <2%; >1 discordant calls in 69 study duplicates; >1 Mendelian error in 9 HapMap trios; Hardy–Weinberg equilibrium P value <10−4 in either white or black subjects; and sex differences in rates of allele frequency ≥0.2 or of heterozygosity >0.3 (see Supplementary Table 2 for SNP filters, available at http://links.lww.com/PAIN/A688). The SNP filters removed 1.4% of SNPs for quality purposes. An additional 5.9% of SNPs that were monomorphic were likewise removed. The final MAF filtering was based on MAF and sample size considerations. The final cleaned data set used for the analysis and imputation included 1.83 million SNPs, or 71.3% of the total. Imputation quality was assessed using the masked/imputed concordance feature of IMPUTE2; all chromosomes exceeded the recommended 95% overall concordance rate. After applying the OEvar filter for MAF and imputation quality, 9.70 million SNPs were retained for the analysis.

2.6. Genetic association analysis

Genome-wide association was tested using PLINK 1.07 software,[32] on discrete genotype and (for imputed SNPs) allele dosage data sets. Logistic regression models were constructed using additive SNP effects and included sex, age, race, and enrollment site. To account for population stratification, genetic ancestry was modeled using the first 3 principal components of the genotypes. The threshold for genome-wide significance was P < 5.0 × 10−8. In addition to the analysis in the full cohort, we also performed tests of association stratified by sex, as TMD has much greater prevalence in females than in males and is believed to have some sex-specific pathophysiologic mechanisms.

To assess the possible effects of population stratification on association results, we examined all associated SNPs using 2 additional approaches. We first subdivided the overall cohort and male and female groups separately into homogeneous racial groups representing European and African-American ancestry. We assessed association effect size and direction in these sex- and race-specific strata using logistic regression, as described previously. We combined the race-specific results by fixed-effects meta-analysis using PLINK. Statistical significance was also assessed within strata and in the meta-analysis, but was considered of secondary importance given the study was not powered to fully explore stratum-specific associations. In a separate approach, we used the software package GEMMA[52] to apply a univariate logistic mixed model analysis to the overall cohort or to sex-stratified subgroups, as appropriate for the observed associations. This model accounted for population stratification by introducing a term describing an n-vector of random effects, using an n × n matrix where each row and column represents each study participant, and each value describes the relatedness between the 2 corresponding individuals. This relatedness matrix was calculated in GEMMA using the “-gk 2” function on a standardized genotype matrix, as standardization is considered to improve performance when SNPs with lower MAF tend to have larger effects. The model also included covariates as described for the logistic regression analyses; significance was assessed by the Wald test.

Power calculations were performed according to the method of Kuo and Zaykin,[17] which estimates the number of true positives expected among the top ranked P values (see Supplementary Table 3 and its legend for further details, available at http://links.lww.com/PAIN/A688). Assuming a truly etiological SNP to have an average odds ratio (OR) = 1.2 with MAF = 0.15, and the number of such true SNPs to be 25, there is 80% chance that one or more of these true positives will score among the top 41 SNPs ranked by order of association P value. For 50 such true loci, there is 80% power that at least one will be found in the top 10 ranked P values. Assumptions about the number of loci and the effect size were guided by distributions for other complex diseases.

2.7. Estimation of heritability

Heritability of susceptibility to TMD was estimated using genome-wide SNP coverage using the GREML method implemented in the software package Genome-Wide Complex Trait Analysis (GCTA).[48] This method estimates the contribution of additive SNP effects to observed variance to estimate narrow-sense heritability or the proportion of phenotypic variance attributable to additive allelic factors. This analysis was performed on a homogenous group of subjects with European ancestry (n = 1627), as indicated by the principal component analysis of the genotypes.

2.8. Meta-analysis

2.8.1. Studies

Seven independent cohorts were used for replication analysis (Table 1). Further details on study design, recruitment, subject characteristics, and phenotyping for each study are provided in Supplementary Materials (available at http://links.lww.com/PAIN/A688).

Table 1

Demographic data for the discovery and replication studies.

Analyses on all replication cohorts were performed with approval from the respective institutional review boards at each site. DNA collected in each replication study was genotyped separately before imputation to the 1000 Genomes Project phase 3 reference panel. Due to differences in genotyping platforms and racial/ethnic composition, quality control procedures were performed separately for each cohort. Tests for association between orofacial pain phenotypes and SNPs were performed using logistic regression, with study-specific covariates. Filtering for MAF and imputation quality was performed using the same OEvar method used for OPPERA.

2.8.2. Meta-analysis statistical methods

Fixed-effects meta-analysis was performed using the metafor package in R version 3.2.2.[45] Allelic log ORs, SEs, and P values were calculated separately for each SNP-phenotype association from each study included in the meta-analysis. The threshold of statistical significance for replication was P < 0.05, divided by the number of independent loci carried forward (ie, 3) using the Bonferroni correction for multiple testing, for the 1-tailed test of association. Study-specific SNP filters for imputation quality and MAF were imposed, as described for each study. The minor allele on the forward strand was used as the effect allele or reference to standardize direction of ORs. The Cochrane's Q statistic P value was used as a measure of heterogeneity of associations among studies.

2.8.3. Functional analysis: expression quantitative trait locus

As ganglia neurons relay sensory information to the brain, they represent a highly relevant tissue in which to study cellular mechanisms of chronic pain. Given the strong correlation in gene expression patterns between trigeminal and dorsal root ganglia (DRG),[3] we examined the correlation of TMD-associated SNPs with gene expression in primary afferent cells from DRG tissue. Dorsal root ganglia specimens were collected postmortem from donors who died following asystole by the CORE/CORID repository at the University of Pittsburgh. Dorsal root ganglia tissues were harvested from bilateral lumbar L4 and L5 and snap-frozen for storage. Thawed samples were genotyped using Illumina's Human Omni Express chip (∼1M probes), while RNA levels were measured with Affymetrix's Human Transcriptome Array 2.0 (∼70K gene-level probes and ∼900K isoform-level probes), which provides for both gene- and exon-level assessment of expression quantitative trait loci (eQTLs) (see diatchenko.lab.mcgill.ca:/DRG-eQTLs/).[28] Both analyses were performed at the Genomics Research Core facility of the University of Pittsburgh. Associations were designated as cis-acting if the distance between the SNP and the gene was less than 1 Mb for gene-level probes, and if the distance was less than 5 kb for exon-level probes. A total of 214 samples with matched genotype/RNA levels were used in this study. We performed eQTL analyses using all samples and then stratified analyses by sex. Donors' median age was 52 years; 199 were white, 105 females and 109 males.

To relate eQTLs observed in DRG specimens (for which TMD status is not available) to RNA regulation in TMD cases and controls, we performed the eQTL analysis using whole blood collected in the OPPERA study. Peripheral blood samples were stored in PAXgene tubes and frozen for subsequent processing. Genotyping was performed on the Illumina Human Omni Express chip (∼1M probes), while RNA levels were measured with Affymetrix's Human Transcriptome Array 2.0 (∼70K gene-level probes and ∼900K isoform-level probes) at the McGill University and Génome Québec Innovation Centre. A total of 167 samples with matched genotype/RNA levels were used. Subjects' median age was 27 years. One hundred seventeen were white, 29 African American, and 9 Hispanic; 115 subjects were females and 52 males. Analyses were performed on the full set of samples as well as in strata subdivided by sex. Genotype differences in probe means were assessed using linear regression to generate P values, and a false discovery rate (FDR) of 1% was set to control for multiple testing.

Although cis-regulation is presumed to result from the direct action of regulatory elements on nearby genes, trans-regulation likely involves downstream effects within biological pathways. We selected expressed genes (eGenes) at FDR 1% level (n = 1231) from P values of association between genes and cis- (nearby genes) and trans- (distal genes) acting eQTLs of interest in human DRG. This data set was submitted to pathway analysis approaches, using QIAGEN's Ingenuity Pathway Analysis (IPA; QIAGEN, Redwood City, CA) and Pathway Studio.[26] Core analysis was applied as a primary IPA analytic tool for identifying cellular pathways, biotoxicology functions, diseases, and known biomarkers associated with the trans-acting gene set. Pathway Studio pathway analyses were performed using the online web server. Analyses were restricted to Gene Ontology (GO) Consortium's biological processes. Statistical significance was assessed based on a FDR of 1% to account for multiple testing and redundancy in gene composition between pathways.

2.8.4. Mras null mutant mouse behavioral assays

A rodent gene-knockout experiment was designed to test functional consequences of genetic variants implicated with TMD in the preceding GWAS and eQTL studies. Mras (muscle RAS oncogene homolog) null mutant mice (B6; 129S5-MrasGt(OST108223)Lex/Ieg) were generated by the Wellcome Trust Knockout Mouse Resource, and frozen embryos were obtained from the European Mouse Mutant Archive (EMMA) repository (EM:02322). Cryorecovery was performed by the Goodman Cancer Research Center Transgenic Core Facility (McGill University), and heterozygous mice were mated in house to generate homozygous, heterozygous, and wild-type mice. All animals used in this experiment were bred on site, and littermates representing the 3 Mras genotypes were introduced into the study as soon as they matured to 8 weeks of age, accounting for the unbalanced numbers between males and females. All mice were housed in standard polycarbonate cages in groups of 3 or 4 same-sex, same-genotype mice in a temperature-controlled (20 ± 1°C) environment (14:10-hour light/dark cycle; lights on at 07:00 hours); tap water and food (Harlan Teklad 8604; Envigo, Huntingdon, United Kingdom) were available ad libitum. Experimenters were blinded to genotype throughout the experiment. All procedures were approved by the Downtown Animal Care Committee at McGill University and conformed to ethical guidelines of the Canadian Council on Animal Care.

All experiments took place during the light cycle, no earlier than 09:00 hours and no later than 16:00 hours. Naive mice were placed in custom-made Plexiglas cubicles (5 × 8.5 × 6 cm) on a perforated metal floor and were habituated for at least 1 hour before testing. The up-down method of Dixon was used to estimate 50% withdrawal thresholds using calibrated von Frey nylon monofilaments[4] (Stoelting Touch Test). Filaments were applied to the plantar surface of the hind paw for 3 seconds and responses were recorded. At least 2 consecutive measures were taken on each hind paw at each time point and averaged. Although data were collected on both hind paws, only data from the hind paw injected with complete Freund's adjuvant (CFA, see below) are presented, as no significant genotype effects on the contralateral paw were observed.

All mice received 20 μL unilateral injections of 50% CFA (Sigma-Aldrich, St. Louis, MO) into the plantar surface of the left hind paw on day 0. Mice were tested for mechanical sensitivity using von Frey fibers before at baseline (1 day before CFA injection) and 3, 7, 10, 14, and 21 days after injection to quantify mechanical allodynia.

Statistical analyses were conducted using an α level of 0.05. Data were analyzed in Systat v. 13 by repeated-measures (RM) 2-way analysis of variance (ANOVA) followed by a post hoc Sidak test with reported P values adjusted for multiple comparisons and graphed with GraphPad Prism v.7.

3. Results

3.1. Discovery genome-wide association scan

After quality control exclusions and genome-wide imputation (see Supplementary Tables 1 and 2, available at http://links.lww.com/PAIN/A688), 999 painful TMD cases and 2031 pain-free controls were analyzed for association with ∼9.7M SNPs (Supplementary Fig. 1, available at http://links.lww.com/PAIN/A688), using logistic regression under an additive model in both the full cohort (Fig. 2A) and in females (Fig. 2B) and males (Fig. 2C) separately. The GREML method implemented in the software package GCTA was used to estimate the variance explained by autosomal SNPs to estimate the narrow-sense heritability of chronic TMD. This analysis indicated that about 17% of the phenotypic variance was accounted for by additive genetic factors.

Figure 2.

Manhattan plots for the test of association with TMD. All imputed SNPs passing the effective sample size filter are included. Threshold lines indicate suggestive (yellow, P = 1 × 10−5) and genome-wide significant (red, P = 5 × 10−8) values; statistically significant SNPs are labeled. (A) All OPPERA subjects, n = 3030. (B) Females only, n = 1956. (C) Males only, n = 1074. SNP, single-nucleotide polymorphism; TMD, temporomandibular disorder.

All SNPs that exceeded the whole genome–corrected threshold for significance are noted in Table 2, including the closest gene features. An imputed insertion/deletion (indel) polymorphism on chromosome 4, rs5862730, achieved genome-wide statistical significance (OR = 1.4, 95% confidence interval [CI]: 1.26-1.61, P = 2.82 × 10−8, Fig. 3A). In a females-only analysis (Fig. 2B), rs5862730 on chromosome 4 was likewise associated with TMD (OR = 1.54, 95% CI: 1.33-1.79, P = 1.7 × 10−8). Another SNP, rs10092633 on chromosome 8, was also significant in females (OR = 4.12, 95% CI: 2.50-6.79, P = 2.9 × 10−8, MAF = 0.04); this SNP appears to be specific to African Americans (minor allele A frequency in whites = 0.0, MAF in African Americans = 0.15; see Supplementary Table 4 for an investigation of potential population stratification effects, available at http://links.lww.com/PAIN/A688).

Table 2

OPPERA TMD genome-wide significant associations.

Figure 3.

Regional plots for genome-wide–associated SNPs. Plots were generated using LocusZoom. The index SNP for each plot is shown as a purple diamond. Linkage disequilibrium (using values obtained from hg19/1000 Genomes EUR data) with other SNPs is indicated by dot color, according to r2 values as indicated in the legend; gray values indicate no pairwise LD information is available for the index SNP. (A) Chromosome 4. The most strongly associated SNP, rs5862730, was not included in the reference panel used for pairwise LD information with neighboring SNPs, so rs59742968 (r2 = 0.95 with rs5862730) was used as the index SNP for this plot. (B) Chromosome 8, females only. (C) Chromosome 3, males only. LD, linkage disequilibrium; SNP, single-nucleotide polymorphism.

In a males-only analysis (Fig. 2C), 3 SNPs on chromosome 3 were significantly associated with TMD (Fig. 3C). These SNPs, which were in high linkage disequilibrium (r2 ≥ 0.8), were rs34612513 (OR = 3.00, 95% CI: 2.06-4.38, P = 1.5 × 10−8); rs28865059 (MAF = 0.09, OR = 2.94, 95% CI: 2.02-4.27, P = 2.2 × 10−8); and rs13078961 (OR = 2.94, 95% CI: 2.02-4.27, P = 2.2 × 10−8). Further statistical analysis revealed a significant interaction between rs13078961 (selected to be illustrative of this locus) and sex (P = 0.0001), and no association signal in the corresponding females-only analysis (OR = 1.02, P = 0.88). Altogether, the data suggest this locus is associated with TMD-only in males.

3.2. Replication and meta-analysis

Five SNPs from 3 distinct loci, significantly associated in either the combined or sex-stratified analyses (Table 2), were carried forward for examination in 7 replication cohorts (Fig. 4). Single-nucleotide polymorphism rs5862730 (chromosome 4) significant in the analysis of all OPPERA subjects was not in the cleaned data sets for SHIP nor NFBC, although it was not significantly associated with TMD at P < 0.017 (1-tailed test, corrected for replication testing at 3 independent loci) in combined or female-only analyses in any of the available replication cohorts (Fig. 4A). A fixed-effects meta-analysis summarizing effects from available replication cohorts (but not the discovery cohort) also did not indicate an overall association effect of this SNP (OR = 1.0, 95% CI: 0.93-1.07).

Figure 4.

Forest plots for associated SNPs in OPPERA and the replication cohorts. The OR and 95% CI for the OPPERA discovery cohort are given in the first row. Odds ratio and 95% CI values for all available replication cohorts are provided in subsequent rows, and the final row indicates the meta-analysis effect in replication cohorts only, not including the discovery cohort. The P value at the top of each plot pertains to the meta-analysis effect in replication cohorts only. (A) chr4: rs5862730, male and female combined; for the replication meta-analysis, Cochrane's Q P = 0.37, I2 = 5.4. (B) chr8: rs10092633, females only, Cochrane's Q P = 0.16, I2 = 38.3. (C) chr3: rs13078961, males only, Cochrane's Q P = 0.79, I2 = 0.0. (D) chr3: rs34612513, males only, Cochrane's Q P = 0.95, I2 = 0.0. (E) chr3: rs28865059, males only, Cochrane's Q P = 0.65, I2 = 0.0. CI, confidence interval; OR, odds ratio; SNP, single-nucleotide polymorphism.

The SNP on chromosome 8, rs10092633, associated with TMD-only in females in the discovery phase, was genotyped in 4 replication data sets (excluding SHIP and NFBC). No significant association with TMD was observed (Fig. 4B), although a suggestive effect in a direction consistent with OPPERA was observed in UKB females (OR = 1.87, 95% CI: 0.83-4.24, 1-tailed P = 7.8 × 10−2), albeit with a low MAF (0.003).

The 3 SNPs on chromosome 3 with male-specific associations were not examined in the SPB cohort, which included only females. Genotypes for 2 SNPs were available in all 6 replication data sets. One SNP, rs13078961, showed nominal replication of the OPPERA association in a meta-analysis of all 6 replication cohorts' males (OR = 1.16, 95% CI: 1.0-1.35, 1-sided test P = 2.3 × 10−2; Fig. 4C), although this result was not significant after Bonferroni correction for the number of independent loci. Among the individual replication cohorts, it was also nominally significantly associated (in males only) in HCHS/SOL (OR = 1.23, 95% CI: 0.98-1.5, 1-sided test P = 3.7 × 10−2; Fig. 4C). The other SNP, rs34612513, had a weaker association in the replication cohorts and did not replicate in the meta-analysis (OR = 1.09, 95% CI: 0.94-1.27, 1-sided P = 1.4 × 10−1; Fig. 4D). Genotypes for one SNP, rs28865059, were present in only 4 of the 6 replication data sets. Similar to rs13078961, the OPPERA finding replicated in the males-only meta-analysis (OR = 1.17, 95% CI: 1.0-1.37, 1-sided test P = 2.4 × 10−2; Fig. 4E) as well as HCHS/SOL males (OR = 1.23, 95% CI: 0.98-1.55, 1-tailed P = 3.7 × 10−2; Fig. 4E), but not in the other 3 cohorts' males. The FDR for replication testing of the 3 discoveries is 0.069, which indicates an expectation of 0.2 false discoveries among the 3 loci.

Thus, only the chromosome 3 SNPs were able to be independently replicated. As replication was nominal and was not significant after Bonferroni correction for the number of tested loci, we further examined this locus for potential functional relevance to pain in the bioinformatic analysis.

3.3. Bioinformatic analyses

To provide evidence for the functional significance of genetic variants that were associated with painful TMD in discovery or replication association studies, we performed further bioinformatics analysis. As disease-associated loci may alter gene function through effects on genes distal from the polymorphic variants identified in the association analysis, we first examined correlations between SNPs and gene expression. Separate analyses considered regulatory effects on expression in cis (ie, with neighboring genes) and globally in trans (with genes located more than 1M nucleotides away); these analyses also considered probes that characterize exons constituting distinct isoforms of genes. The same 5 SNPs that emerged from the discovery cohort were investigated for eQTLs in both human postmortem nerve tissue extracted from DRG and in blood from OPPERA participants.

In the DRG experiment, none of the 5 SNPs were associated in cis with mRNA expression at the gene level, using a FDR threshold of 1%. However, at the exon level, the 3 SNPs associated with TMD in males-only (rs34612513, rs13078961, and rs28865059) exhibited cis-acting eQTL effects (Supplementary Table 5, available at http://links.lww.com/PAIN/A688). These SNPs showed significant association with expression of the gene muscle RAS oncogene homolog (MRAS, representative SNP rs13078961: beta = −0.51, P = 2.43 × 10−5, FDR = 5.83 × 10−3; Fig. 5), which has a transcription start site about 380 kb from the eQTL region (Fig. 3C). In sex-stratified analyses, significant associations between rs13078961 and several MRAS exons were observed in males (for exon 8: beta = −0.55, P = 2.25 × 10−4, FDR = 4.73 × 10−3) but not in females (beta = −0.47, P = 8.2 × 10−3, FDR = 2.0 × 10−1).

Figure 5.

Replicated male-only TMD risk SNP rs13078961 is an eQTL for the 3′UTR of MRAS gene. (A) MRAS has 6 different splice isoforms in the RefSeq database (top), comprising 9 exons (numbered 1-9 inclusively). mRNA levels measured by 9 nonoverlapping microarray probes in exons 5 to 9 were significantly associated with dosage of minor alleles of SNP rs13078961 in human DRG as discovery (green), blood as replication (blue), and TMD caseness (red). (B) The relative position of microarray probes at 3′ UTR of MRAS (exon 9) associated with SNP rs13078961 in human DRG (green), blood (blue), and TMD caseness (red). Coordinates are UCSC hg19. For brevity, only the last 2 digits of probe identifiers are shown; whole probe identifiers are: PSR030135xx. DRG, dorsal root ganglia; eQTL, expression quantitative trait locus; SNP, single-nucleotide polymorphism; TMD, temporomandibular disorder; UTR, untranslated region.

For trans-acting associations, a large number of genes showed statistically relevant associations at FDR 1% with the same male-specific locus. This association was observed at both gene and exon levels, and in both cases was driven exclusively by males (Supplementary Table 5c, available at http://links.lww.com/PAIN/A688). To assess whether identified trans-acting genes are enriched in any biological pathway or are related to known human diseases, we applied a pathway analysis approach. Both Ingenuity Pathway Analysis and Pathway Studio methods identified immune response and particularly B Cell and T Cell Receptor Signaling as of highest significance (Supplementary Tables 6 and 7, available at http://links.lww.com/PAIN/A688).

Next, we examined eQTL association using blood RNA samples from the OPPERA cohort. Similarly to the DRG eQTL analysis, the representative SNP rs13078961 was associated with MRAS mRNA expression at the exon level with a similar pattern in males and females combined (beta = −0.21, P = 3.76 × 10−2), and nominally significant only in males in a sex-stratified analysis (beta = −0.39, P = 2.17 × 10−2) (Supplementary Table 8, available at http://links.lww.com/PAIN/A688). Importantly, in all cases, the minor allele of rs13078961 was associated with a lower level of MRAS expression. We found no statistically significant associations at FDR 1% levels for trans-acting eQTLs for either gene- or exon-level probes. Thus, our data indicate MRAS mRNA exhibits a similar cis-eQTL–dependent expression pattern in sensory DRG and whole blood in males.

Finally, we tested whether whole blood mRNA levels of MRAS varied according to TMD case status. We found a statistically significant association at the gene level (beta = −0.04, P = 4.8 × 10−2) in males and females combined. One exon of MRAS in males was also nominally significantly associated with a much larger effect size (beta = −0.27, P = 1.7 × 10−2) (Supplementary Table 9, available at http://links.lww.com/PAIN/A688). The direction of association was consistently negative, indicating levels of mRNA are lower in TMD cases, in agreement with the eQTL results.

Overall, the male-specific eQTL and differential M-RAS expression was observed mostly with different exons of the MRAS transcript, rather than with the entire gene (Fig. 5). These probes, although distinct, are situated adjacent to each other at the very end of the 3′ untranslated region (UTR), which is the most common location of eQTL regulation of mRNAs.[28,33] Importantly, we note that we did not apply Bonferroni correction for the number of exons we tested, but such correction would be overly conservative, as expression of the vast majority of exons is not independent, and notably, our association pattern was consistently located at the last 3′ UTR exon of the gene.

3.4. Mras null mutant mouse nociceptive assays

Genetic association and eQTL analyses both suggested MRAS expression level influenced chronic pain in males; so, we tested a null mutant (knockout) mouse in which functional Mras protein was reduced or eliminated. The injection of CFA was selected as an assay relevant to a wide variety of chronic pain conditions, especially those with an inflammatory neuroimmune component, to confirm the involvement of this gene in pain. We compared the mechanical hypersensitivity (allodynia) induced by hind paw injection of CFA between Mras null mutant, heterozygote, and wild-type mice (Fig. 6). All groups exhibited a decrease in withdrawal threshold to stimulation with von Frey fibers by day 3, which resolved by day 21 (Fig. 6A). There were no genotype differences in baseline mechanical thresholds before CFA injection or peak mechanical allodynia observed 3 days after CFA injection. A 2-between (genotype × sex) RM ANOVA (baseline to day 21) showed a significant genotype × sex × time interaction (F10,120 = 2.1, P = 0.03). In sex-segregated one-between RM ANOVAs, there was a significant genotype × RM interaction in males (F10,80 = 2.9, P = 0.004; Fig. 6B) but not in females (F10,40 = 2.1, P = 0.06; Fig. 6C). In male but not female null mutant mice, the duration of CFA-induced mechanical allodynia was significantly prolonged, lasting for 14 to 21 days instead of the 3- to 7-day duration exhibited by wildtypes.

Figure 6.

Mechanical allodynia after CFA injury in Mras knockout mice and wild-type controls. Mean ± SD reported. (A) All mice (n = 30); (B) male mice only (n = 19); and (C) female mice only (n = 11). *P < 0.05; **P < 0.01. CFA, complete Freund's adjuvant.

4. Discussion

Extending beyond previous candidate gene investigations, we now report a discovery GWAS and validation study of clinically diagnosed painful TMD. In the discovery phase analysis of 3030 individuals from the OPPERA study, genome-wide significant associations with chronic TMD were found for 3 loci, including an association in males-only at 3 linked SNPs on chromosome 3 (rs28865059, rs13078961, and rs34612513), with OR around 3.0. In meta-analysis of up to 7 independent case–control and population-based cohorts, we observed a nominal replication of 2 SNPs in this locus. Functional studies of the chromosome 3 locus implicated allele-specific expression differences of the nearby gene MRAS for its effects, and we observed that male Mras knockout mice exhibited slower recovery from inflammatory pain hypersensitivity induced by CFA.

As in many chronic idiopathic pain conditions,[23] a higher prevalence of TMD among females has been observed, suggesting that sex-specific mechanisms may be responsible. We examined genotype effects within sex to discover variants and pathways potentially underlying these mechanisms. Relevant associations at the chromosome 3 locus were observed in males-only, and the sex-specificity of this effect was supported by a significant sex × genotype interaction (P = 0.0001). Consistent with the primary association of these SNPs with TMD in males-only, all eQTL associations were nominally significant in males in both DRG and blood RNA samples. Sex-dependent effects observed in the mouse inflammatory pain model provide insight into the nature of the susceptibility to persistent pain conferred by a decrease in Mras expression. Neither sex showed a genotype difference in baseline mechanical pain sensitivity, and all mice displayed similar levels of acute pain hypersensitivity 3 days after CFA treatment. In males, null mutant mice showed substantially delayed recovery from CFA at day 10 and 14 compared with wildtypes, suggesting that low expression of Mras interferes with recuperative processes. Although the mechanism by which Mras facilitates this recovery is unknown, our results indicate that it is an important component to healing in males, whereas females appear to rely on separate pathways.

Our results demonstrate the importance of conducting sex-stratified analyses in association and behavioral studies of pain.[24,25,44] We identified more genome-wide significant associations when we separated our cohorts into males and females, compared with the combined cohort, despite reduced sample size.[2,46] Our findings suggest that males and females may develop TMD through different pathophysiological pathways. An interesting question remains as to why the strongest association results, which we were able to reproduce and show functional relevance, are in male-specific loci, although chronic TMD is much more prevalent in females. In contrast to the typical approach of explaining increased risk in females, perhaps these findings highlight the value of understanding lowered risk or resiliency in males for pain chronicity. In this regard, our discovery may have identified a pathway that protects males from developing chronic TMD.

M-RAS is a member of the Ras family of small GTPases and is broadly expressed in many tissue types including brain, nerves, and heart. It is generally involved in cell growth and differentiation processes, and is a component of tumor necrosis factor-alpha (TNF-α)[49] and ERK-MAPK[50] signaling pathways. Because the candidate SNPs were associated with the last exon's expression level rather than the entire MRAS mRNA, we suggest that these SNPs or haplotypes affect exon expression in a relatively less abundant alternative isoform, by creating or removing alternative exons, the polyadenylation signal, or splice sites,[18,22,51] or via selective miRNA-dependent degradation.[20]

The direction of the eQTL effect was identical for all 3 associated SNPs at this locus, such that the minor allele associated with higher risk of TMD was also associated with lower MRAS expression in all exons in both DRG and blood (Supplementary Tables 5 and 8, available at http://links.lww.com/PAIN/A688). Consistent with this observation, we found mRNA levels of MRAS negatively correlated with TMD case status with the strongest effect again in the probe SR03013530 at the 3′ UTR (Supplementary Table 9, available at http://links.lww.com/PAIN/A688; Fig. 6), reinforcing the direct contribution of MRAS to the TMD-resistant phenotype.

As MRAS has not previously been connected with a pain phenotype or disorder, we sought to generate further hypotheses regarding the downstream effects of the MRAS eQTL using pathway analysis methods. Our analysis of trans-acting associations with the 3 associated SNPs further confirmed that the primary link of this locus with TMD was exclusive to males, as multiple trans-associated genes were found in our male-specific DRG analysis. We used 2 separate pathway analysis programs to understand whether these genes are related to particular biological processes and identified immune response with the focus on B/T-cell receptor signaling as a primary pathway in DRG. This finding further emphasizes a role of the immune system in neurological tissues in the pathophysiology of chronic pain and is in line with our recent evidence showing contribution of the human leukocyte antigen locus to pain-related genetic associations.[28] It is also consistent with the growing understanding of the sexually divergent roles of immune cells, as animal studies have provided evidence that pain hypersensitivity in males is dependent on microglia, whereas females rely on a separate system likely involving T lymphocytes.[42,43] Sex differences in M-RAS signaling may arise due to this dimorphism in male and female systems, or they may be secondary to interactions between TNF-α and hormones such as estrogen[40] or testosterone.[5,36] Given that we identified much stronger SNP-dependent MRAS regulation in DRG than in blood (Supplementary Table 5, available at http://links.lww.com/PAIN/A688), our pathway analysis results suggest that this association reflects interaction of the immune system with sensory neurons (ie, neuroimmune modulation), rather than immune processes per se.

Other findings from the discovery study implicated genes with the potential for intriguing pathophysiological insights, but until they can be reproduced in independent cohorts should be considered provisional. We were not able to identify any relevant functionality for rs5862730, the indel on chromosome 4, through our bioinformatics approaches, nor were we able to conclusively pinpoint the gene responsible for the effect. The most promising nearby candidate gene is SMAD1, which is a transcriptional modulator involved in cell growth and development, most prominently through mediation of signaling by bone morphogenetic proteins.[16] Although a clear link between SMAD1 and TMD pain has not been confirmed, other SMAD pathway genes have been linked to acute inflammatory (but not neuropathic) pain[31] and excitability of DRG neurons[14] as well as osteoarthritis of the temporomandibular joint.[47] Furthermore, our group found that TGFβ1 regulation of cytokines such as MCP-1 and IL-8, which is mediated by SMAD protein activation, is consistently inhibited in the plasma of TMD cases relative to controls.[39]

As with all genome-wide surveys, this study has a number of limitations. An estimated 88% of GWAS findings are located outside of gene-coding regions,[9] a proportion consistent with the findings of this study. We have used bioinformatics tools to determine functionality of our findings, leading to the implication of the putative eQTL at the “males-only” chromosome 3 locus with MRAS expression. However, the mechanisms underlying the remainder of the putative associations must be elucidated by further study. This study was also marked by the low replication rate typically observed in GWAS, commonly ascribed to low power, as well as phenotypic or ancestral heterogeneity.[12,15] We have attempted to replicate statistically significant findings from a moderately powered discovery study (n = 3030) in cohorts of vastly larger size; however, greater sample size often comes at a huge cost to phenotypic homogeneity. Although we accounted for genetic ancestral variation within studies through principal component analysis adjustment, regional and ethnic differences between the replication cohorts may alter genetic effects, precluding replication. Several of the SNPs observed in the discovery cohort had relatively low MAF in particular racial groups, which would limit our power to observe effects in other cohorts, especially those with different ancestral background. Although we observed a number of significant associations, the associated SNPs account for a very small proportion of the observed phenotypic variance. As we estimated the additive heritability of TMD at 17%, there are likely many TMD susceptibility variants yet to be discovered with larger cohorts and improved phenotyping. Finally, the consequences of Mras deletion in the knockout mouse model may be quite different from the allelic effect on regulation of MRAS in humans, limiting the interpretability of the animal study.

As has been observed in GWAS meta-analyses of other complex traits, we anticipate studies combining larger cohorts in diverse populations will increase the number of loci associated with TMD. Our results provide evidence to guide the design of these studies, suggesting that different molecular mechanisms underlie pathology in males and females. These results also provide targets for investigation into novel therapeutic strategies, notably the Ras subfamily of signaling molecules, and further emphasize the contribution of immune processes in the pathophysiology of chronic pain.

Conflict of interest statement

L. Diatchenko and W. Maixner are cofounders and equity stock holders in Algynomics, Inc, a company providing research services in personalized pain medication and diagnostics. S.B. Smith, G.D. Slade, and R.B. Fillingim are equity stock holders in Algynomics, Inc. The remaining authors have no conflicts of interest to disclose.