Methods
MCP was found to be a complex, polygenic trait genetically correlated with psychiatric and other disorders in a 2019 GWAS ( 5 ). Recent changes to ICD-11 coding for chronic pain and International Association for the Study of Pain definitions of chronic pain ( 33 – 35 ) support the study of chronic pain as a disease. Genes involved in central nervous system and immune function were found to be associated with MCP using MAGMA ( 36 ), and gene expression of MCP-related genes was enriched in the brain. Summary statistics were used for transcriptome-wide association study analysis through the TI approach S-PrediXcan ( 13 ).
GREX was imputed using MCP GWAS output and TI models from the GTEx (Genotype-Tissue Expression Project) ( 37 ) in 13 brain tissues ( Table S1 ) using S-PrediXcan. Multiple testing correction (Bonferroni) was applied and resulted in 2 thresholds for significance: 1) a per-tissue threshold correcting for all genes tested in each tissue ( Table S1 ), and 2) an experiment-wide threshold correcting for all genes across all tissues ( p = 7.9 × 10 −7 ). Then, we sought to replicate our findings using a different TI method, summary–transcriptome-wide association study ( 38 ) (see the Supplement ).
A recent genetic study of pain intensity was carried out in 598,339 Million Veterans Program participants ( 39 ) and included FUSION transcriptome-wide association analysis and prediction models for 6 brain tissues (anterior cingulate cortex, cerebellar hemisphere cortex, frontal cortex, cerebellum, and dorsolateral prefrontal cortex). Pain intensity was significantly genetically correlated with MCP ( r g = 0.79) ( 39 ). We downloaded the 361 significant gene-tissue results [Supplementary Table 20 in Toikumo et al. ( 39 )] and carried out a Fisher’s exact test to ascertain whether results overlap represented significant replication.
Pathway analyses were carried out using FUMA GENE2FUNC ( 40 ) including all per-tissue significant gene results ( n = 89). We tested for enrichment of all gene sets available in FUMA GENE2FUNC with all genes that had at least one S-PrediXcan prediction model available and were included in FUMA as background ( n = 15,588). Significant gene results were also investigated using the FUMA DrugBank (see the Supplement ).
We queried Connectivity Map (CMap), a large database of perturbation signatures maintained by the Broad Institute ( 41 , 42 ), using genes up- and downregulated in MCP ( Table S2 ). We filtered results to retain compounds (drugs) passing CMap quality control with significant connectivity scores (−log 10 (FDR [false discovery rate]–corrected p ) > 1.3, FDR–corrected p < .05).
To probe relationships between MCP-associated GREX and clinical phenotypes, we performed a series of PheWASs (see the Supplement ) in the Mount Sinai Bio Me biobank.
Bio Me is a large, diverse, hospital-based biobank that includes electronic health record and genotype data for 31,704 participants in the first data freeze. A total of 1236 phecodes for Bio Me participants were included in analyses presented in this paper. Phecodes are a high-throughput method that reduce electronic health record dimension and complexity in which ICD-10 codes are manually grouped according to clinical similarity ( 43 ). Here, we used previously curated phecodes ( 44 ). A full list of phecodes can be searched at https://phewascatalog.org/phecodes_icd10 or through download of the “PheCode Definitions v1.2 ICD-10-CM map” available at https://phewascatalog.org/phecodes_icd10cm .
First, we imputed MCP-GREX (chronic pain–related genetically regulated gene expression) for 31,704 Bio Me freeze 1 participants, split across 6 genotype-derived ancestry groups ( Table S3 ).
Specifically, we imputed GREX in all 13 brain regions and in whole blood for all 89 unique genes previously identified as significant MCP-GREX. We tested for associations between these GREX values and Bio Me phecodes with at least 10 available cases in at least one ancestry [total phecodes = 1236 ( 44 )]. Results were meta-analyzed using inverse variance-weighted meta-analysis in METAL ( 45 ). Multiple testing correction (within-gene FDR) was then applied.
To validate our MCP associations, we tested whether MCP-associated genes were associated with pain. A numeric rating scale (NRS) ranging from 0 to 10, where 10 is the worst pain possible and 0 is no pain, was recorded for Bio Me participants and aggregated into a mean pain score across instances in which the pain NRS was recorded. Associations were tested between significant MCP-GREX results and mean pain scores, and results were meta-analyzed across ancestry groups using inverse variance-weighted meta-analysis in METAL. FDR correction was performed as previously described.
Results
We applied S-PrediXcan to the largest available summary statistics for MCP ( n = 387,649). We identified 95 experiment-wide significant gene-tissue associations ( p < 7.97 × 10 −7 ), including 36 unique genes ( Table 1 ). An experiment-wide threshold is likely overly conservative because many expression quantitative trait loci are shared between tissues; therefore, we also applied a within-tissue Bonferroni threshold ( Table S1 ; Figure 1A , B ). We identified an additional 134 gene-tissue associations that reached within-tissue significance, including 53 additional unique genes.
Of these 89 genes, 59 were not previously associated with MCP ( 5 ) ( Table S4 ; Figure S2 ). We also found significant levels of replication of our gene-tissue findings in summary–transcriptome-wide association study ( Supplement; Tables S5 , S6 ). We also found significant replication of S-PrediXcan findings within significant TI findings for pain intensity. Six significant gene-tissue associations for MCP ( Tables S7 , S8 ) were also significant in analyses of pain intensity, representing significant replication ( p = 4 × 10 −9 ). To test whether significant associations were enriched in specific brain regions, we compared the proportion of experiment-wide significant associations per region with the proportion of genes tested in that region (binomial enrichment tests). We found significantly more experiment-wide significant associations in the nucleus accumbens basal ganglia than would be expected by chance (14.7% vs. 7.6%, p Binomial = .0075) and significantly fewer in the cerebellar hemisphere (4.2% vs. 9.0%, p Binomial = .038). Repeating this test for nominally associated genes, 3 brain regions showed fewer associations than would be expected by chance: the hippocampus (5.3% vs. 5.8%, p Binomial = .033), spinal cord cervical C1 (4.4% vs. 5.1%, p Binomial = .0014), and substantia nigra (3.4% vs. 4.0%, p Binomial = .0035).
To identify functional patterns of MCP-GREX associations, we conducted a gene set enrichment analysis using FUMA (see the Supplement ). Genes associated with MCP-GREX were significantly enriched in the positional gene set chr3p21 ( p = 5.27 × 10 −19 ) ( Figure 2A ), which was also implicated in anorexia nervosa ( 46 ). MCP-GREX genes were significantly enriched for genes associated with 8 GWASs ( Figure 2B ). This included a previous GWAS of MCP ( p = 5.54 × 10 −6 ) ( 5 ), sleep duration (short sleep) ( p = 2.27 × 10 −11 ), extremely high intelligence ( p = 6.66 × 10 −8 ), regular attendance at gyms and sports clubs ( p = 6.66 × 10 −8 ), and religious group attendance ( p = 7.66 × 10 −6 ), as well as inflammatory conditions (ulcerative colitis, p = 1.95 × 10 −5 , inflammatory bowel disease, p = 5.9 × 10 −3 ) and age at first birth ( p = 1.57 × 10 −3 ). FUMA DrugBank lookups ( Table S9 ) identified 19 genes as drug targets. CMap analyses identified 23 compounds with significant connectivity scores ( Table 2 ).
To probe clinical consequences of our MCP-associated genes, we performed a PheWAS in the Mount Sinai Bio Me biobank. First, we imputed MCP-GREX for 89 significant MCP-GREX gene-tissue associations for 18,806 biobank participants who had available mean pain score data and tested for association between GREX and mean pain score. We identified 37 associations including 10 unique genes between MCP-GREX and mean pain score ( Table 3 ). Next, we tested for phenome-wide associations, imputing MCP-GREX for 89 significant MCP-GREX gene-tissue associations for 31,704 Bio Me participants across 6 ancestry groups. Then, we meta-analyzed across ancestry using METAL and applied multiple testing correction (FDR). We identified 16 significant GREX-phecode associations across 9 brain regions, including 10 unique gene-phecode associations ( Table 3 ; Figure 3 ). Associated phecodes included cardiac dysrhythmia, metabolic syndrome, disc disorders/dorsopathies, joint/ligament sprain, anemias, and neurological disorders.
Because pain and chronic pain are core symptoms of many of these diagnoses, and some genes with significant MCP-GREX were significantly associated with pain NRS, it is difficult to discern whether our MCP genes are associated with pain experience or directly with the trait itself. Therefore, we repeated our PheWAS on a subset of Bio Me participants and included mean pain scores derived from pain NRS information as covariates. We also carried out a PheWAS with adjustment identical to our main analyses (no adjustment for mean pain score) on the same subset of participants. We found the results to be significantly different from the main PheWAS results, but after comparison with the unadjusted-subset PheWAS, this appears to have been driven by a reduction in sample size rather than by mean pain score ( Tables S10 , S11 ). Sample size is significantly reduced when adjusting for pain score because many Bio Me participants do not have pain NRS information available.
Discussion
These results reveal novel genes, theoretically enriched for causal effect, that are relevant to chronic pain development, thus providing new insight into mechanisms of chronic pain. By applying TI using S-PrediXcan, we were able to perform a well-powered study of gene expression in brain tissue and whole blood, which is currently not feasible with existing cohorts in which chronic pain phenotyping, genotype, and expression data are available together due to limited sample sizes. In the following section, we contextualize our findings with a focus on MCP-GREX genes found to be significantly associated with clinical traits (phecodes) in our Bio Me PheWAS analysis.
GREX of ILRUN , involved in innate immune response and highly expressed in B cells ( 47 ), was significantly associated with MCP in the basal ganglia of the nucleus accumbens, hypothalamus, amygdala, and cortex in the original S-PrediXcan analysis and with primary thrombocytopenia across all 4 tissues in our PheWAS ( Table 4 ). Primary thrombocytopenia is an autoimmune platelet disorder that causes low peripheral plate counts and symptoms including joint and abdominal pain, bleeding, and bruising. ILRUN has also been linked to the renin-angiotensin-aldosterone system (involved in blood volume, sodium reabsorption, and vascular tone among other processes) in a study of SARS-CoV-2 infection ( 48 ). Peripheral small Ad and C fibers that transmit pain signals contain cells expressing renin-angiotensin-aldosterone system components, and renin-angiotensin-aldosterone system modulators have been shown to affect pain relief ( 49 ). Our results suggest a role for ILRUN in the brain in chronic pain development, in addition to in pain perception in the periphery.
MCP-GREX of MON1B in both the amygdala and cervical spinal cord C1 was found to be significantly associated with anemias ( Table 4 ); this phecode includes sickle cell anemia, thalassemia, and hemolytic anemias, all of which have often been associated with significant pain ( 50 ). Iron deficiency and iron-deficiency anemia are also generally associated with chronic inflammatory disease and chronic pain ( 51 ). Dysregulation of iron metabolism can play a key role in immune cell homeostasis and inflammation ( 52 , 53 ). MON1B also encodes a protein for which defects are associated with autoimmune pathology ( 54 ), a process that plays a significant role in chronic pain ( 55 ). This protein is also a key regulator of endocytic sorting by Numb, and so is linked to cell migration, asymmetrical cell division, and differentiation ( 56 ).
DCAKD encodes a protein linked to neurodevelopment ( 57 ) that is expressed widely in the brain ( 58 ), and MCP-GREX of this gene in the caudate basal ganglia was negatively associated with cardiac dysrhythmia ( Table 4 ). Previous studies indicate a relationship between magnetic resonance imaging markers of cerebral small vessel disease and DCAKD ( 59 ) and Friedrich’s ataxia ( 60 ), a disease of progressive neurodegeneration, heart, and spinal problems ( 61 , 62 ). Heart rate variability is thought to represent hyperarousal and has been linked to emotion regulation and chronic pain ( 63 , 64 ). In addition, certain nerve blocks can treat both cardiac and chronic pain conditions ( 65 ).
ECM1 encodes a protein involved in type 2 helper T cell migration ( 66 ) and skin development ( 67 ). In PheWAS analyses, ECM1 MCP-GREX was associated with dysmetabolic syndrome X (aka metabolic syndrome) in 3 different brain tissues ( 68 , 69 ) ( Table 4 ). This syndrome has been associated with increased risk of cardiovascular disease and type 2 diabetes ( 68 , 70 ). T cells have been associated with insulin resistance development in obesity ( 71 ); having metabolic syndrome can affect T cell development [reviewed by ( 72 )]; and the amount of memory T cells has been associated with a proinflammatory state ( 73 ). These cell types could be therapeutic targets in chronic pain treatment ( 74 – 77 ) and could represent a sex-dimorphic mediator of pain hypersensitivity [reviewed by ( 78 , 79 )].
PACSIN3 encodes a protein involved in the actin cytoskeleton and formation of vesicles ( 80 ). This protein also binds TRPV4; channelopathy mutations in the TRPV4 gene lead to skeletal dysplasias, Charcot-Marie-Tooth disease subtype 2C, premature osteoarthritis, and neurological disorders ( 81 ). TRPV4 channels are also important in skin function ( 82 ) and are involved in the itch-scratch cycle ( 83 , 84 ). TRP channels have also been implicated in chronic low back pain ( 85 ) and investigated as a therapeutic target in fibromyalgia ( 86 , 87 ). PACSIN3 MCP-GREX in the basal ganglia of the nucleus accumbens was significantly associated with bullous dermatoses in PheWAS analyses ( Table 4 ). Bullous dermatoses are autoimmune skin conditions of painful blistering ( 88 – 91 ). Although itch and pain are considered to be distinct ( 84 ), they share many similarities ( 92 ). Results here suggest that TRPV4 ion channels and their interaction with PACSIN3 could be a point of overlap between chronic pain and itch.
RAD51 is involved in DNA repair ( 93 , 94 ). RAD51 mutations have been linked to congenital mirror movement disorder ( 95 ) and cancers ( 96 ). MCP-GREX at this gene in the substantia nigra was significantly associated with disturbances of sulfur-bearing amino acid metabolism ( Table 4 ). This phecode includes homocystinuria (the body is unable to process methionine) and methylenetetrahydrofolate reductase (MTHFR) deficiency (homocysteine levels are elevated) ( 97 ). Both processes are part of DNA metabolism ( 98 ), and elevated homocysteine levels are associated with a range of illnesses and neurotoxicity ( 99 ). RAD51 foci (indicators of cellular replication stress) ( 100 ) were increased in experiments examining folate deficiency ( 101 ). Previous studies in rodents showed that elevated homocysteine caused mechanical allodynia ( 102 ), and PheWAS results indicate a role for this mechanism of sensitization in human chronic pain.
SCAI encodes a transcriptional cofactor that regulates invasive cell migration ( 103 ), including in gliomas ( 104 ). MCP-GREX of this gene in the cortex was associated with toxic/inflammatory neuropathy in PheWAS analyses, and this gene was differentially expressed in rat models of diabetic neuropathy in the spinal cord ( 105 ). Our findings suggest a similar role for human SCAI in neuropathy.
SLC38A3 encodes a glutamine transporter ( 106 ) involved in cell energy metabolism. Glutamine is the preferred energy source for rapidly proliferating cell populations in the nervous system, immune system, and cancer cells ( 107 – 111 ). SLC38A3 is also expressed in muscles, and significant MCP-GREX in the caudate basal ganglia was found to be associated with joint and muscle sprain ( Table 4 ), suggesting that the glutamine transporter encoded by SLC38A3 has a central as well as a peripheral role. SLC38A3 MCP-GREX in the same brain area was also significantly associated with neurological disorders ( Table 4 ), consistent with research showing relationships between glutamine metabolism in the brain and neurological conditions ( 112 – 115 ). GABAergic (gamma-aminobutyric acidergic) gene regulatory elements have also been implicated in neurological and psychiatric diseases ( 116 – 120 ), glutamate receptors in neurological dysfunction ( 121 ), and treating neurodegeneration through targeting glutamate transporters ( 122 ). Activity-dependent synaptic plasticity also involves glutamate and glutamine metabolism ( 123 ). Glutamine has also been investigated as a chronic pain biomarker because concentrations vary in individuals with chronic pain compared with control participants ( 124 , 125 ), and glutamine supplementation may be helpful in vaso-occlusive crisis in sickle cell disease ( 126 ). Glutamine levels have also been associated with individual pain sensitivity differences ( 127 ) and migraine ( 128 ). Finally, glutamine supplementation was associated with reduced opioid use in sickle cell disease in a small study, highlighting potential as a harm- and pain-reducing compound in chronic pain treatment ( 129 ). Finally, ERICH2 MCP-GREX in the amygdala was significantly associated with dorsopathies ( Table 4 ).
Psychiatric disorder–related phecodes and phecodes assigned to chronic pain conditions, e.g., rheumatoid arthritis or endometriosis, were not significantly associated with MCP-GREX. In contrast, significant genetic correlations between MCP and, e.g., major depressive disorder and MCP and rheumatoid arthritis were found in a previous study ( 5 ). Genetic correlations are calculated using all single nucleotide polymorphism associations genome wide rather than at a gene level, which may explain these differences. In addition, in theory, S-PrediXcan results represent gene expression changes that occur before chronic pain development (whereas GWAS summary statistics used in linkage disequilibrium score regression represent genetic associations more generally). This suggests that the gene expression changes that contribute to chronic pain development do not directly contribute to psychiatric conditions (e.g., major depressive disorder), which is consistent with previous studies that have suggested that chronic pain can have a causal effect on major depression development but not vice versa ( 5 ). Another possibility is that tissues that were not examined in this study are associated with MCP-GREX and would show associations with psychiatric disorder or other expected phecodes in a PheWAS. However, it is difficult to explain why these nonbrain tissues, and not brain tissue, would show this result. We chose to examine brain and whole blood because chronic pain involves significant changes in the brain and spinal cord ( 16 – 19 ), and whole blood represents a tissue of interest due to immune components and ease of testing for, e.g., potential chronic pain biomarkers. Finally, phecodes generally represent a broad category of diagnoses; for example, the phecode for mood disorder (296) encompasses depression associated with major depressive disorder, bipolar disorder, and schizophrenia, and this heterogeneity could affect PheWAS results.
After adjusting our PheWAS association testing for mean pain score, results were significantly different compared with the main PheWAS analyses. However, these changes appear to be driven by reduction in sample size because unadjusted and adjusted analyses in the same subset of individuals showed similar results. Although NRS is a widely used pain reporting measure in clinical and research settings ( 130 ), it can change in unpredictable ways over time in chronic pain ( 131 , 132 ), may not accurately reflect treatment outcome when used alone ( 133 ), and may not be the most useful measure for identifying clinically important pain ( 134 ) or changes in pain ( 135 ). Pain NRS may not represent an ideal assessment tool in nonacute pain at the population or group level despite some studies demonstrating stability when an NRS was used to assess improvement in individuals over time ( 136 ) because perception of pain, which influences NRS ratings, is likely to be significantly different between individuals with and without chronic pain ( 137 ). People with chronic pain may rate moderate to high levels of pain as tolerable ( 138 ); conversely, depression or depressive symptoms that are commonly comorbid with chronic pain could lead to the reporting of higher NRS scores ( 139 – 141 ).
Chronic pain is complex and difficult to treat successfully. The results shown here could inform treatment development; genes where MCP-GREX is associated with upregulation may present better targets in genomic medicine (downregulation of a gene can be easier to induce than upregulation), and genes where significant MCP-GREX is shown in a singular tissue may present a better target for potential animal modeling of chronic pain compared with genes where MCP-GREX is widespread. DrugBank lookups provide suggestions for drug repurposing, and several drugs highlighted are already used experimentally in chronic pain treatment, e.g., monoclonal antibodies in migraine ( 142 – 144 ) and drugs that increase inhibitory glycinergic neurotransmission in the spinal cord ( 145 , 146 ). Several compounds identified in CMap analysis also show potential in chronic pain treatment; PX-12 showed anti-allodynia effects in a rodent model of chronic pain ( 147 ); physostigmine showed an antihyperalgesic effect in clinical trials ( 148 ); and SR-2640 activates TREK-1 channels that are associated with nociceptive hypersensitivity in rodent models ( 149 ). Arcyriaflavin-a is a potential therapeutic compound in endometriosis ( 150 ), as sorbinil ( 151 ) and fenoterol ( 152 ) are in diabetic neuropathy. Ursolic acid has demonstrated antinociceptive properties in animal models ( 153 ), and analgesic properties of palmitoylethanolamide ( 154 ) and luteolin ( 155 ) have been shown in multiple studies. Other findings are established pain treatments, e.g., aspirin and nimesulide. Other compounds, e.g., epidermal growth factor receptor (EGFR) inhibitor PD-153035, affect cancer-related pathways, which are also implicated in chronic pain ( 156 ), thus presenting novel treatment targets.
Conclusions
We carried out the largest TI study of a chronic pain trait to date, making important progress in translating GWAS findings into insights into chronic pain development and beginning to bridge the gap between genotype (GWAS output) and phenotype (MCP). Specific brain tissues and the direction of effect of MCP-GREX are also given; pathways of interest and potential mechanistic overlap with other medical conditions are indicated; and several genes showing significant MCP-GREX are also potential drug targets. We also identified several compounds with opposite expression perturbation signatures to MCP (i.e., potentially therapeutic compounds in chronic pain). Results of our PheWAS in which we adjusted for mean pain score indicate that associations tend not to be driven solely by pain perception. PheWAS results indicate potential shared causal pathways between chronic pain and conditions such as metabolic syndrome, anemias, and cardiac dysrhythmia.