A multi-ancestry genetic study of pain intensity in 598,339 veterans

Abstract

Chronic pain is a common problem, with more than one-fifth of adult Americans reporting pain daily or on most days. It adversely affects the quality of life and imposes substantial personal and economic costs. Efforts to treat chronic pain using opioids had a central role in precipitating the opioid crisis. Despite an estimated heritability of 25-50%, the genetic architecture of chronic pain is not well-characterized, in part because studies have largely been limited to samples of European ancestry. To help address this knowledge gap, we conducted a cross-ancestry meta-analysis of pain intensity in 598,339 participants in the Million Veteran Program, which identified 126 independent genetic loci, 69 of which are new. Pain intensity was genetically correlated with other pain phenotypes, level of substance use and substance use disorders, other psychiatric traits, education level and cognitive traits. Integration of the genome-wide association studies findings with functional genomics data shows enrichment for putatively causal genes (n = 142) and proteins (n = 14) expressed in brain tissues, specifically in GABAergic neurons. Drug repurposing analysis identified anticonvulsants, β-blockers and calcium-channel blockers, among other drug groups, as having potential analgesic effects. Our results provide insights into key molecular contributors to the experience of pain and highlight attractive drug targets.

Extended Data Fig. 1 |. Overview of the study.

Top left: primary GWAS analyses for pain intensity. Within ancestry, GWAS for African American (AA), European American (EA) and Hispanic American (HA) followed by cross-ancestry meta-analysis. These results were used for all downstream analyses. Top right: secondary GWAS analyses for pain intensity. Bottom: downstream analyses were conducted using the cross-ancestry, AA and EA GWAS results as indicated by color shadings: primary GWAS (green) and supplementary GWAS (brown).

Extended Data Fig. 2 |. Manhattan plot for the pain intensity in European American GWAS analysis.

Identified 87 independent risk loci. Novel loci (n = 52) are annotated in pink. The red line indicates GWS after correction for multiple testing (P < 5 × 10⁻⁸).

Extended Data Fig. 3 |. Effect-effect plot of cross-ancestry meta-analyses lead SNPs in the primary and secondary GWASs.

The magnitude and direction of the effect sizes are plotted for each GWAS. The results show significant (P < 2.2 × 10⁻¹⁶) high correlation (Pearson r test, two-sided) between the effect sizes (β) of pain intensity lead SNPs for primary GWAS and those for non-OUD (r = 1, a), non-zero (r = 0.97, b), males (r = 1, c) and females (r = 0.88, d).

Extended Data Fig. 4 |. LDSC genetic correlations for pain intensity primary and secondary GWAS.

African American: primary GWAS, n = 112,968; non-OUD GWAS, n = 104,050; non-zero GWAS, n = 61,499; male GWAS, n = 97,343; female GWAS, n = 15,625. European American: primary GWAS, n = 436,683; non-OUD GWAS, n = 416,740; non-zero GWAS, n = 202,784; male GWAS, n = 404,510; female GWAS, N = 32,173. Error bar is presented as 95% confidence interval.

Extended Data Fig. 5 |. MAGMA tissue enrichment for pain intensity in cross-ancestry and European American GWAS results.

Tissue enrichment analyses were conducted using FUMA. Bonferroni correction threshold (represented by the black dashed line) = 9.25 × 10⁻⁴ (0.05/54).

Extended Data Fig. 6 |. Gene-based Manhattan plots for cross-ancestry, European American and African American GWAS.

Gene-based association analyses were conducted using FUMA and genes that survive multiple correction are annotated (Bonferroni p = 2.67 × 10⁻⁶ [0.05/18,702]).

Extended Data Fig. 7 |. Regional plot for TRAIP*rs2247036 and MST1R*rs9815930 on chromosome 3.

Credible locus prioritized by FINEMAP (PP > 0.5) is annotated with red rings. The MST1R*rs9815930 locus is in high LD (r² > 0.8) with the lead variant TRAIP*rs2247036.

Extended Data Fig. 8 |. Regional plot for NOP14*rs71597204 and GRK4*rs2798303 on chromosome 4.

Credible locus prioritized by FINEMAP (PP > 0.5) is annotated with red rings. The GRK4*rs2798303 locus is in moderate LD (r² > 0.4) with the lead variant NOP14*rs71597204.

Fig. 1 |. Manhattan plot for the pain intensity cross-ancestry GWAS meta-analysis (n = 598,339).

This identified 125 independent index variants. The nearest gene to the 66 new loci (65 autosomal and 1 X-chromosomal) is annotated. SNPs above the red line are GWS after correction for multiple testing (P < 5 × 10⁻⁸).

Fig. 2 |. Enrichment of pain intensity in the brain.

a, Partitioning heritability enrichment analyses using LDSC showing enrichment for pain intensity in the CNS, adrenal, liver, cardiovascular, skeletal and immune/hematopoietic tissues. The dashed black lines indicate Bonferroni-corrected significance for multiple testing (P < 0.005). b, Proportion of heritability shows robust enrichment for SNPs in brain and immune-related tissues. c–e, Heritability enrichment analyses for gene expression (c and d) and chromatin interaction (top 35 annotations are shown in e; see Supplementary Table 17 for full details) using Genotype-Tissue Expression (GTEx) data show enrichment for pain intensity in brain regions previously associated with chronic pain. Bonferroni correction was applied within each tissue conditioned on the number of genes tested.

Fig. 3 |. Gene prioritization for pain intensity.

a, Genomic annotation of credible sets using FINEMAP shows enrichment largely in noncoding regions and to a lesser extent in exons. PIP, posterior inclusion probability; UTR, untranslated region. b, Annotation of known and new credible genes. Dashed line indicates PP > 0.5. c, Number of overlapping genes across functional prediction models. d, Tissue enrichment of prioritized genes using SMR and GTEx data shows enrichment in brain regions. The size of the circle reflects −log₁₀(P). Bonferroni correction was applied within each tissue conditioned on the number of genes tested.

Fig. 4 |. Genetic correlation.

Genetic correlation for pain intensity using LDSC. All points passing Bonferroni correction (P = 5.62 × 10⁻⁴ (0.05/89)) are plotted. The color of the circle indicates the phenotypic category. The vertical dashed line represents genetic correlation = 0. BMI, body mass index; HDL, high-density lipoprotein; T2D, type 2 diabetes; CUD, cannabis use disorder; CWP, chronic widespread pain; PAU, problematic alcohol use; MDD, major depressive disorder; ADHD, attention deficit/hypersensitivity disorder.

Fig. 5 |. Drug repurposing.

Druggable targets and drug interactions for eight credible genes associated with pain intensity. For a full list of credible drug targets, see Supplementary Table 37.