Study overview

Genome-wide discovery meta-analysis was comprised of adults of European ancestry from 16 cohorts (n = 158,025 including 29,531 CBP cases; Table 1), including 15 CHARGE cohorts and participants from the UK Biobank (UKB) interim data release (UKB1). After quality control, the number of SNPs included in the meta-analysis ranged from 6,205,227 to 9,775,703, depending on the cohort (S1 and S2 Tables). Linkage disequilibrium (LD) score regression (LDsr) was used distinguish polygenicity from potential confounding,[18] using LD scores from European ancestry 1000 Genomes data. The genome-wide significance level was defined as p<5×10⁻⁸, and suggestive significance level was defined as p<5×10⁻⁷, after using the LDsr intercept as a correction factor. For those SNPs of genome-wide suggestive significance in the discovery phase, replication was conducted in a sample of UKB European ancestry participants (UKB2) who were not part of the interim data release (n = 283,752 subjects, including 50,915 CBP cases), and a joint (discovery-replication) meta-analysis was performed. We then conducted functional characterization of variants and loci achieving genome-wide significance in the joint meta-analysis.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Cohorts in meta-analysis of genome-wide association studies of chronic back pain (European ancestry).

https://doi.org/10.1371/journal.pgen.1007601.t001

Meta-analysis of GWAS of CBP

The characteristics of cohorts included in the discovery meta-analysis are shown in Table 1. The mean age of participants in each cohort ranged between 53–76 years. Within cohorts, mean age, BMI, and proportion of females was more often higher among CBP cases than among those without CBP. A quantile-quantile plot comparing the meta-analysis association results with those expected by chance is presented in S1 Fig. The LDsr intercept was 1.007 (standard error 0.006), λ was 1.114, and the LDsr ratio was 0.0581 (standard error 0.053), providing no evidence of inflation of p-values from population stratification. Meta-analysis results are summarized in the Manhattan plot shown in Fig 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Manhattan plot for meta-analysis (discovery) of GWAS of chronic back pain (n = 158,025).

GWAS = genome-wide association study. Results use the linkage disequilibrium score regression (LDSR) intercept as a correction factor. Red line depicts genome-wide statistical significance (P <5×10⁻⁸). Blue line depicts suggestive significance (P <5×10⁻⁷).

https://doi.org/10.1371/journal.pgen.1007601.g001

A genome-wide significant association (OR 1.08, p = 7.2×10⁻¹⁰) was found for rs12310519 on chromosome 12 in an intronic region of SOX5, with little evidence for heterogeneity (I² = 0, p = 0.95) (Table 2, S2 Fig). Several other signals were in high LD (r²>0.8) with the top signal (S3 Fig), but none were independently associated with CBP in analyses conditional on rs12310519.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Association results for chronic back pain: Meta-analysis (discovery), replication, and joint meta-analysis^*.

https://doi.org/10.1371/journal.pgen.1007601.t002

No other variants achieved genome-wide significance, but variants in three other loci reached suggestive significance (Table 2, S3 Table, S4–S9 Figs): rs1453867 (OR 0.95, p = 7.7×10⁻⁸), located in an intronic region of chromosome 2 within DIS3L2; rs7833174 (OR 1.06, p = 1.0×10⁻⁷), located in an intergenic region on chromosome 8 between CCDC26 (a long non-coding RNA) and GSDMC; and rs4384683 (OR 0.95, p = 3.2×10⁻⁷), located in an intronic region of chromosome 18 within DCC. In each of these 3 regions, there was no other variant reaching the suggestive significance level in analyses conditional on the lead SNP in the region. Post hoc secondary analyses of the discovery sample showed effects of similar magnitude and direction between the CHARGE cohorts and the UKB interim data release for associations between the lead variants in the top 4 loci and CBP (S4 Table).

We examined these 4 top variants in 283,752 UKB individuals not included in the discovery sample (UKB2), including 50,915 cases (Table 2). For all 4 variants, the direction of association was the same in discovery and replication. The association for rs12310519 in SOX5 replicated in UKB2 (OR 1.06, p = 5.3×10⁻¹¹), and exceeded genome-wide significance in the joint analysis (OR = 1.07, p = 4.5×10⁻¹⁹). Of the 3 suggestive-significance variants from the discovery stage, rs7833174 at CCDC26/GSDMC (OR 1.05, p = 4.4×10⁻¹³) and rs4384683 in DCC (OR 0.97, p = 2.4×10⁻¹⁰) exceeded genome-wide significance in the joint meta-analysis, but rs1453867 in DISL32 (OR 0.98, p = 3.9×10⁻⁷) did not (Table 2). Thus, we demonstrate genome-wide significant associations of CBP with loci tagged by rs12310519 (SOX5), rs7833174 (CCDC26/GSDMC), and rs4384683 (DCC), with replication for rs12310519 in SOX5.

Characterization of variants in SOX5, CCDC26/GSDMC, and DCC

Functional characterization followed the same steps for each of the 3 loci that achieved genome-wide significance in the joint meta-analysis. First, we examined cross-phenotype genetic associations between each lead SNP and traits with possible conceptual links to CBP, in look-ups of publicly and privately available GWAS datasets. Where the lead SNP was not present in a dataset, we examined associations with the variant in highest LD with the lead SNP. Second, we annotated lead variants and those in LD (r²≥0.6) for consequences on gene functions (using the combined annotation dependent depletion [CADD] score [21]), potential regulatory functions (using RegulomeDB score[22]), and effects on gene expression (using GTExv6 [23, 24]), and examined whether these variants resided in enhancer regions for selected tissues with connections to the spine or pain processing (using data from the Roadmap Epigenomics Consortium [25, 26]) (Methods, S1 Text).

SOX5.

Among CBP-related traits examined, the lead SNP in SOX5, rs12310519, was most strongly associated with imaging-detected lumbar intervertebral disc degeneration (p = 1.1×10⁻⁴; S5 Table)[27]. The highest CADD score among SOX5 variants was 10.52 for the lead SNP in the region rs12310519, indicating it is predicted to be among the 10% most deleterious possible substitutions in the human genome (S1 Appendix). However, the overall regulatory potential of these variants was low according to RegulomeDB score (scores of 6 [‘minimal binding evidence’]) (S1 Appendix). There were no meaningful associations with gene expression using GTExv6. The lead SNP rs12310519 and variants in LD (r²>0.6) contained active enhancer marks in chondrogenic cells using Roadmap Epigenomics Consortium data (Fig 2, S1 Appendix).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Significant variants are co-localized with potential gene regulatory markers.

The heatmap depicts the percentage of variants in gene regulatory regions (associated with enhancers/promotors) in LD (r²>0.6) with rs7833174 (CCDC26/GSDMC), rs12310519 (SOX5), and rs4384683 (DCC). We examined epigenetic histone marks in selected cell types/tissues including chondrogenic, bone-related, neuronal and brain cells/tissues; mesodermal cells (related to notochord); and psoas muscle (located proximal to the lumbar spine).

https://doi.org/10.1371/journal.pgen.1007601.g002

Secondary analyses to examine interrelationships with height

Post hoc analyses among UKB1 participants from the discovery sample (n = 120,023) indicated that associations with CBP for the lead variant in SOX5 were similar with and without adjustment for height as a covariate, and in conditional analyses accounting for height (S1 Text, S7 Table). Associations with CBP for the lead variant in CCDC26/GSDMC were also similar with and without adjustment for height as a covariate. However, associations with CBP were markedly diminished when conditional on the lead height-associated variant in the region, and associations with height were markedly diminished when conditioned on the lead CBP-associated variant in the region (S7 Table). This suggests that the same functional variant is responsible for association of CCDC26/GSDMC locus with height and CBP, although an alternative explanation is two functional variants in tight LD. Associations with CBP for the lead variant in DCC were similar with and without adjustment for height as a covariate, and in conditional analyses accounting for height (S7 Table).

To examine possible causal effects of height on CBP, we conducted a two-sample Mendelian randomization (MR) analysis using genetic variants associated with standing height in the GIANT consortium as the exposure, and the discovery phase GWAS meta-analysis of CBP as the outcome. Results of the two-sample MR using 326 SNPs and the instruments involved are available in the S2 Appendix. These instruments explained 10.1% of the variance in height, with an average SNP-height F-statistic of 78.5, indicating substantial instrument strength. ORs for CBP were 1.10 per standard deviation increase in height with inverse variance weighted (IVW) regression (p<0.0001). However, there was significant heterogeneity among SNPs (I² = 0.35; p<0.0001), suggesting horizontal pleiotropy for at least some SNPs (S3 Appendix). Estimates with other MR methods were directionally consistent, and all but MR-Egger regression were statistically significant: an OR for CBP of 1.09 per standard deviation increase in height with MR-Egger regression (p = 0.19); 1.14 per standard deviation increase in height with the weighted median method (p<0.0001), and 1.17 per standard deviation increase in height with the weighted mode method (p = 0.02) (S3 Appendix). The magnitude of MR estimates after excluding 14 outlier SNPs were very similar to the two-sample MR using 326 SNPs, but with substantially less heterogeneity amongst SNPs (I² = 0.12; p = 0.04), just exceeding nominal significance (S2 and S3 Appendices). MR-Egger intercepts were close to 0 with both the 326 SNP and 312 SNP instruments, and neither were statistically significant, suggesting no strong directional horizontal pleiotropy under the InSIDE (Instrument Strength Independent on Direct Effect) assumption (S3 Appendix).

Heritability of CBP and genetic correlations

SNP heritability of CBP on the liability scale was 7.6%. Partitioned heritability by functional category using stratified LD score regression showed significant enrichment (p = 0.0004) for regions conserved in mammals, with 2.6% of SNPs explaining 40% of the SNP heritability of CBP, without significant enrichment for other functional categories, including coding regions (S4 Appendix). This pattern of partitioned heritability was broadly similar across cell type groups, including central nervous system (CNS), connective tissue and bone, and skeletal muscle, among others (S4 Appendix). Genetic correlations of nominal significance (range of r_g 0.17–0.31, p<0.05) were found with anthropometric traits involving obesity or body fat distribution (waist circumference, hip circumference, waist-hip ratios, overweight/obesity classes, and BMI), but not with height (S9 Table). Larger magnitude nominally significant (p<0.05) genetic correlations were also found with depression-related phenotypes (range of r_g 0.46–0.52), self-reported osteoarthritis (r_g = 0.63), and ICD-10-defined osteoarthritis (r_g = 0.49) phenotypes.

Study design and populations

Discovery meta-analysis included adults of European ancestry from 16 population- and community-based cohorts: Cardiovascular Health Study (CHS), Framingham Heart Study (FHS), Generation Scotland (GS), Johnston County Osteoarthritis Project (JoCo), Osteoporotic Fractures in Men (MrOS) Sweden (MrOS-Gothenburg and MrOS-Malmo), MrOS US, Osteoarthritis Initiative (OAI), Rotterdam Study (RS1, RS2, and RS3), Study of Osteoporotic Fractures (SOF), 10,001 Dalmatians (Vis and Korcula), TwinsUK, and UKB participants from the interim data release[20]. Replication was conducted among UKB European ancestry participants not included in the discovery stage (UKB2), and a joint (discovery-replication) meta-analysis was performed. The separation of analyzed data from UKB into discovery (UKB1) and replication phases (UKB2) reflects the history of this scientific collaboration, in that our initial meta-analysis plan included only the UKB data then available to us and for which we had obtained approval to use (UKB1 [the interim data release]). By the time our meta-analysis was completed, all UKB data had become available; the remainder of UKB data was therefore used for replication. Detailed descriptions of the study cohorts are provided in Table 1 and the Supplemental Methods. This meta-analysis was approved by the Research and Development Committee of VA Puget Sound Health Care System (RDIS 0010, MIRB 00903). Institutional Review Board/Ethics Committee approvals at the individual study sites include those listed in the S2 Text. Written or electronic consent was provided for all studies.

Chronic back pain (CBP) phenotype

There is no “gold-standard” definition for CBP. Consistent with the most commonly accepted clinical and research definitions for CBP [3, 4], CBP cases were defined in this study using one of 3 definitions depending on the cohort (Table 1, S10 Table): 1) ≥3 months of back pain, 2) ≥6 months of back pain, and 3) ≥1 month of back pain in consecutive years (reflecting ≥12 months of back pain). For each cohort, the comparison group (“controls”) was comprised of those who reported not having back pain or reported back pain of insufficient duration to be included as a case. This study used a general definition examining chronic ‘back pain’, as opposed to a more specific chronic ‘low back pain’ definition, due to the fact that most of the included cohorts did not include question items permitting localization of pain to the low back or lumbar region specifically.

Genotyping

Details of genotyping, quality control, imputation methods, and genome-wide analysis for each cohort were study-specific (S1 and S2 Tables). In brief, genotyping was performed using commercially available genome-wide arrays. Imputation of single nucleotide polymorphisms (SNPs) and insertions/deletions (indels) was performed using reference panels from 1000 Genomes phase 1 version 3 or phase 3,[71] or the Haplotype Reference Consortium.[72] Analyses of UKB participants was restricted to the White British ancestry subset who self-report as White, British, and have very similar genetic ancestry backgrounds based on the results of principal components analysis (PCA); further quality control followed recommended practices for UKB[73] (S1 Table).

Statistical analysis

We conducted genome-wide association analyses in each of the 16 cohorts, and subsequent meta-analysis of autosomal SNPs to combine results from all cohorts. Each site conducted GWAS using logistic regression models with additive genetic effects to test for associations between each variant and CBP as a binary trait. These models adjusted for age, sex, study-specific covariates, and population substructure using principal components (S2 Table). Height and body mass index (BMI), calculated as weight in kilograms divided by height in meters squared, were not included as covariates in site-specific GWAS, since these traits might lie along the causal pathway or (in the case of BMI) reflect a consequence of CBP. Harmonization and quality control of GWAS results from each cohort were conducted using the EasyQC software package in the R statistical environment (v3.2.2), using methods described previously.[74] After removal of SNPs with low minor allele frequencies (<0.005 for UKB, <0.03 for Vis, <0.01 for other cohorts) or imputation quality (<0.7 for UKB, <0.6 for other cohorts), deviation from Hardy-Weinberg equilibrium (p < 1 x 10–6), low number of cases (<15) or controls (<15), large absolute values of beta coefficient (≥10), and low minor allele count (≤10), call rate <0.95, the range of SNPs included in the meta-analysis was between 6,205,227 (Croatia-Vis) and 9,775,703 (MrOs-Gothenburg) (S2 Table). Fixed-effect inverse-variance weighted meta-analysis was performed with METAL version 2011-03-25 (http://csg.sph.umich.edu/abecasis/metal/), using the LDsr intercept as a correction factor. The meta-analysis was filtered for variants with fewer than 125,000 informative participants, to ensure that SNP-CBP associations were informed by a plurality of cohorts, and not only the UKB interim data release. For this reason, only variants with MAF<0.01 (SNPs) were included in the meta-analysis. Quality control and meta-analysis were conducted twice, independently of each other, by researchers at the University of Washington (MP and PS) and at PolyOmica (YT, YA, and LCK). The results from the two centers were compared to ensure accuracy. Q-Q and Manhattan plots were generated in R. We conducted conditional and joint (COJO) analysis using summary data (http://cnsgenomics.com/software/gcta/#About) to examine associations conditional on the most significant variant at each locus (S1 Text).

The most highly-associated variants at genome-wide significant loci were subjected to replication among UKB participants not included in the discovery sample (UKB2). Analysis in UKB2 used logistic regression with additive genetic effects, adjusting for age, sex, array, and principal components (S2 Table); significant replication was defined using a Bonferroni-corrected threshold of p<0.05 divided by the number of genome-wide significant loci. The most highly-associated variants at loci with suggestive significance were selected for a joint (discovery-replication) meta-analysis using p<5×10⁻⁸ to define genome-wide significance. Further details regarding analysis are provided in S1 Text. A post hoc analysis was conducted to stratify the discovery phase meta-analysis by the CHARGE cohorts (meta-analysis of 15 GWAS) vs. the UKB interim data release cohort (S4 Table).

For genome-wide significant variants, we examined GWAS associations with selected traits with possible links to CBP (anthropometrics, arthritis, depression and depressive symptoms, and imaging-based spinal degeneration) in publicly and privately available GWAS datasets. We conducted functional annotation using FUMA (http://fuma.ctglab.nl). FUMA draws upon multiple publicly available databases, annotating variants for consequences on gene functions using the combined annotation dependent depletion (CADD) score,[21] potential regulatory functions (RegulomeDB score),[22] and effects on gene expression using expression quantitative trait loci (eQTLs) of different tissue types (GTExv6) [23, 24] (S1 Text). The higher the CADD score, the more potentially deleterious is the variant. A CADD score of ≥10 indicates a variant predicted to be among the 10% most deleterious substitutions involving the human genome, a score of ≥20 indicates a variant among the 1% most deleterious, and so forth.[21] We used data from the Roadmap Epigenomics Project to evaluate whether the lead variants at each locus and those in LD (r²>0.6) reside in enhancer regions for selected tissues with possible conceptual connections to the spine via roles in chondrogenesis, vertebral development, muscle, and pain processing in the CNS.[25, 26].

Because two CBP-associated variants were found to be associated with height in prior published GWAS, we conducted post hoc region-specific secondary GWAS analyses accounting for height, among UKB participants from the discovery stage. Further details of functional annotation and secondary analyses are provided in S1 Text. We also conducted a two-sample Mendelian randomization to examine potential causal effects of height on CBP using significant variants associated with standing height from the GIANT consortium as the exposure, and the discovery phase meta-analysis of CBP, using the R package MRbase. We used the inverse-variance weighted regression (IVW) approach as our primary analysis method,[75] and additional analyses with other MR methods (MR-Egger regression, weighted median function, and weighted mode); presenting the results yielded from different MR methods is recommended to demonstrate sensitivity to different patterns of assumption violations.[63, 76] We examined heterogeneity among SNPs using forest plots, funnel plots, heterogeneity statistics, and the MR-Egger intercept test for directional horizontal pleiotropy. Further details of MR methods are provided in the S1 Text.

Given that 30% of UKB participants are related to at least one other person in the cohort, we also conducted post hoc secondary GWAS in UKB1 to examine whether relatedness might have influenced our results. These analyses used linear mixed-effect models (BOLT-LMM), adjusting for age, sex, study-specific covariates, and principal components. The statistical significance of GWAS results for UKB1 using BOLT-LMM were descriptively compared with the original results using logistic regression, for the lead variants achieving suggestive significance in the GWAS meta-analysis.

Finally, we used LDsr of summary-level GWAS results from the discovery stage to estimate heritability due to common autosomal SNPs and genetic correlations.[18] We transformed the observed SNP heritability to the liability scale, in order to make heritability estimates for CBP comparable with traditional heritability estimates from twin studies.[64] We used stratified LDsr to partition heritability across functional categories of the genome, using methods described previously.[65] The threshold for determining the statistical significance of 53 functional categories in partitioning heritability was set at p<9.4 x 10–4 (0.05/53). Further details of methods for partitioning heritability are provided in the S4 Appendix. We used cross-trait LDsr and publicly available meta-GWAS results from LDhub to examine genetic correlations with selected traits with possible links to CBP: anthropometrics (height, waist/hip circumference, BMI, and overweight/obesity), depression and depressive symptoms, osteoarthritis, and rheumatoid arthritis.[77, 78]