Genome-wide meta-analysis of 158,000 individuals of European ancestry identifies three loci associated with chronic back pain

1Back pain causes more years lived with disability than any other health condition worldwide. Although most adults experience a new (‘acute’) episode of back pain at some point in their lives, the societal burden of back pain is driven by the minority of individuals who fail to recover from such episodes and go on to develop persistent (‘chronic’) back pain. Chronic back pain (CBP) has a number of definitions but is most often considered as back pain of duration ≥3 months in clinical practice, and a duration of ≥6 months is also commonly used in research.

2Back pain is moderately heritable. Meta-analysis of 11 twin studies of back pain indicates a heritability of 40%, with a pattern of monozygotic (rMZ = 0.56) and dizygotic (rDZ = 0.28) twin correlations suggesting an additive genetic model (2rDZ = rMZ). Heritability is greater for chronic than for acute back pain. Nevertheless, genetic studies attempting to identify specific genetic markers for CBP have to date been limited to small studies using the candidate gene approach. Although CBP is often attributed to anatomic changes such as intervertebral disc degeneration or disc herniation, such findings have only weak association with CBP, and explain only a small proportion (7–23%) of the genetic influence on back pain. The unexplained genetic contribution to CBP may involve not only spine pathology but also functional predisposition to chronic pain involving higher-order neurologic processes related to the generation and maintenance of pain. Furthermore, psychological factors such as depression are widely recognized as important risk factors for CBP. Given the range of processes that might contribute to CBP, the agnostic genome-wide association approach may identify novel genetic variants associated with CBP and provide insights into underlying biological mechanisms not previously considered.

3This research was an international collaboration between investigators from the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Consortium Musculoskeletal Workgroup and the European Union FP7 project Pain-OMICS (‘Multi-dimensional omics approach to stratification of patients with low back pain’). We conducted a meta-analysis of GWAS of CBP in adults of European ancestry from 16 community- and population-based cohorts, including those from the CHARGE and PainOmics consortia, and the UK Biobank. The aim was to identify novel associations between specific genetic markers and CBP, and elucidate the biological mechanisms underlying this condition.

4Genome-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, using LD scores from European ancestry 1000 Genomes data. The genome-wide significance level was defined as p<5×10−8, and suggestive significance level was defined as p<5×10−7, 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.

5The 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.

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

7No 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−8), located in an intronic region of chromosome 2 within DIS3L2; rs7833174 (OR 1.06, p = 1.0×10−7), 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−7), 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).

8We 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−11), and exceeded genome-wide significance in the joint analysis (OR = 1.07, p = 4.5×10−19). Of the 3 suggestive-significance variants from the discovery stage, rs7833174 at CCDC26/GSDMC (OR 1.05, p = 4.4×10−13) and rs4384683 in DCC (OR 0.97, p = 2.4×10−10) exceeded genome-wide significance in the joint meta-analysis, but rs1453867 in DISL32 (OR 0.98, p = 3.9×10−7) 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.

9Functional 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 (r2≥0.6) for consequences on gene functions (using the combined annotation dependent depletion [CADD] score), potential regulatory functions (using RegulomeDB score), and effects on gene expression (using GTExv6), 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) (Methods, S1 Text).

10Among 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−4; S5 Table). 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 (r2>0.6) contained active enhancer marks in chondrogenic cells using Roadmap Epigenomics Consortium data (Fig 2, S1 Appendix).

11In look-ups of GWAS of selected traits with possible links to CBP (S5 Table), the lead SNP rs7833174 was most strongly associated with height in UKB (p = 1.3×10−59) and hip circumference in UKB (p = 1.8×10−5). The lead SNP rs7833174 was also associated with radiographic hip osteoarthritis (p = 4.9×10−4) (S5 Table). All variants in CCDC26/GSDMC that were suggestively associated with CBP showed cross-phenotypic associations with lumbar microdiscectomy for sciatica in a recent GWAS of Icelandic adults (S6 Table, lowest p = 5.6×10−12). The direction of effect on these other phenotypes was the same as the direction of effect on CBP (i.e. the T allele associated with greater height, greater risk of radiographic hip osteoarthritis, and greater risk of lumbar discectomy for sciatica was also associated with greater CBP risk). The highest CADD score among CBP-associated variants at CCDC26/GSDMC was 18.75 for rs6470778, indicating that this SNP is predicted to be among the 5% most deleterious substitutions in the human genome, and the overall regulatory potential of these variants was substantial according to Regulome DB score (highest RegulomeDB score of 2b [‘likely to affect binding’]) (S1 Appendix). In examination of effects on gene expression using GTExv6, variants in LD with rs7833174 (r2>0.6) were also eQTLs for GSDMC expression in esophageal mucosa, skin, and skeletal muscle (S1 Appendix, p<5×10−8). The lead SNP rs7833174 and variants in LD (r2>0.6) contained active enhancer marks located in regulatory regions for mesodermal cells, astrocytes, chondrogenic cells, and osteoblasts in data from the Roadmap Epigenomics Consortium (Fig 2, S1 Appendix).

12In look-ups of GWAS of selected traits with possible links to CBP (S5 Table), the lead SNP rs4384683 was associated with depressive symptoms(p = 5.9×10−4), with the same direction of effect (i.e the A allele was associated with less depressive symptoms and lower CBP risk). The highest CADD score among variants at DCC that were suggestively associated with CBP was 11.21 for rs2116378, indicating that rs2116378 is predicted to be among the 10% most deleterious substitutions in the human genome, but the overall regulatory potential of these variants was low according to Regulome DB score (highest RegulomeDB score of 5) (S1 Appendix). There were no meaningful associations with gene expression using GTExv6. In data from the Roadmap Epigenomics Consortium, the lead SNP rs4384683 and variants in LD (r2>0.6) contained active enhancer marks in H9 human embryonic stem cell-derived neural cells (Fig 2, S1 Appendix).

13Post 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).

14To 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 (I2 = 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 (I2 = 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).

15GWAS analyses using linear mixed-effect models in UKB1 yielded associations between the 3 lead SNPs and CBP that were very similar to the original analyses using logistic regression in terms of statistical significance, indicating no meaningful influence on the study results due to relatedness (S8 Table).

16SNP 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 rg 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 rg 0.46–0.52), self-reported osteoarthritis (rg = 0.63), and ICD-10-defined osteoarthritis (rg = 0.49) phenotypes.

17This study is the first meta-GWAS of CBP. This collaboration between two international consortia for genomic studies of complex traits in the USA and Europe incorporated data from 16 cohorts and more than 441,000 participants of European ancestry across discovery and replication samples. Our study identifies three novel associations with CBP for loci at SOX5, CCDC26/GSDMC, and DCC.

18CBP was most strongly associated with rs12310519 in an intronic region of the SOX5 gene. The SOX genes are a family of transcription factors involved in virtually all phases of embryonic development, and are thought to determine the fate of many cell types. The SOX genes are defined by containing the HMG (‘high mobility group’) box of a gene involved in sex determination called SRY (‘sex determining region’). SOX5 and SOX6 have overlapping functions and work together in close coordination that is necessary for efficient chondrogenesis. Inactivation of SOX5 leads to minor defects in cartilage and skeletogenesis in mice, whereas SOX5/SOX6 double knockouts have severe chondrodysplasia. Together with SOX9, SOX5 and SOX6 are sometimes referred to as the ‘master chondrogenic SOX trio’. Prior work indicates an important role for SOX5 in articular cartilage and osteoarthritis, and such a role was also supported by our functional annotation showing that rs12310519 (and SNPs in high LD) overlapped with potential regulatory regions for chondrogenic cells. SOX5/6 are also essential for notochord development, and through this role they are critical in the formation of the vertebral column, including the intervertebral discs. Inactivation of SOX5 and/or SOX6 in mice leads to a range of abnormalities in the development of spinal structures. Although variants in SOX5 have not been reported in prior GWAS of limb osteoarthritis (knee, hip, or hand), the association of SOX5 with CBP may involve the spinal structures specifically. Consistent with this, we found a cross-phenotypic association for the lead CBP-associated variant in SOX5 with imaging-detected lumbar intervertebral disc degeneration in a prior GWAS meta-analysis. Future GWAS may also be useful to characterize other spine-related phenotypes besides disc degeneration, such as osteoarthritis of the zygapophyseal (‘facet’) joint, the only true synovial joint in the spine.

19The intergenic variants at CCDC26/GSDMC associated with CBP in the current study were also previously found associated with lumbar microdiscectomy for sciatica due to intervertebral disc herniation. These findings are intriguing, given that lumbar disc herniations (an aspect of lumbar disc degeneration) have long been implicated as a cause of some forms of back pain. Recent studies have concluded that associations between imaging-detected lumbar disc herniation and CBP are of modest magnitude. This might explain the small magnitude association of the top variant at CCDC26/GSDMC with CBP in the current study (OR 1.08 in discovery), in contrast to the larger magnitude association seen with microdiscectomy for sciatica (OR 1.23). Functional characterization of these intergenic variants suggest the likely involvement of the gene GSDMC. GSDMC encodes the protein Gasdermin C, part of the GSDM family of genes that is expressed in epithelial tissues. Although the specific role of GSDMC in lumbar disc herniation and/or sciatica is unclear, GSDMC is associated with differential methylation patterns in osteoarthritis-related cartilage and subchondral bone cartilage, consistent with our findings that variants in LD with rs4384683 were located in potential regulatory regions in chondrocytes and osteoblasts. Our examination of univariate cross-phenotypic genetic associations for CBP-associated variants at CCDC26/GSDMC also suggest pleiotropy with radiographic hip OA for rs6470763. Taken together, these data suggest interconnections between variants at CCDC26/GSDMC and CBP involving cartilage, osteoarthritis, and/or lumbar disc degeneration.

20The third significant CBP-associated variant in our study was rs4384683, an intronic variant in the gene DCC (Deleted in Colorectal Carcinoma), which co-localized with regulatory regions in neural embryonic stem cells. DCC encodes a transmembrane protein that is a receptor for netrin-1, an axonal guidance molecule involved in the development of spinal and cortical commissural neurons. Interactions between DCC and netrin-1 are among the best-studied axonal guidance processes, with key roles during development and in adulthood, and they also affect angiogenesis. Increased expression of netrin-1 and DCC occurs in degenerate human intervertebral discs compared to healthy control discs, and in nucleus pulposus compared to annulus fibrosis. Netrin-1/DCC might therefore mediate neurovascular ingrowth into the intervertebral disc, which has long been implicated as a possible mechanism of chronic discogenic back pain. Given the well-known phenotypic correlation between depression and CBP, however, another possible explanation for the link between CBP and DCC (suggested by the cross-phenotype association of rs4384683 with depressive symptoms) is pleiotropy. Netrin-1/DCC interactions are also known to play a role in pain processing in the spinal cord in animal models of mechanical allodynia. Taken together, this information suggests various potential mechanisms underlying the association between DCC and CBP, including nociceptive pathways and/or the involvement of mood.

21Some epidemiological studies report that greater height confers increased risk of back pain, although a systematic review found no association. Variants in CCDC26/GSDMC associated with CBP in our meta-analysis were also reported to be associated with height in prior GWAS; hence, post hoc analyses devoted special attention to the role of height in CBP. These region-specific analyses showed that the association of SOX5 and DCC variants with CBP was independent of height; however, CBP- and height-associated variants at CCDC26/GSDMC were tightly linked and could not be disentangled in conditional analyses (S7 Table), indicating that the association of variants in CCDC26/GSDMC with both CBP and height might be explained by biological pleiotropy or mediated pleiotropy (i.e. pleiotropy due to causal effects). Mendelian randomization analysis, drawing on information from hundreds of genetic markers distributed across the genome, suggested that height may have causal effects on CBP, although with a degree of heterogeneity suggesting horizontal pleiotropy for some SNPs. Such evidence of horizontal pleiotropy is common in MR studies of complex traits, and can be seen even in MR studies of exposure-outcome relationships where causal effects are known. Taken together, our findings suggest a causal component to the relationship between height and CBP, but do not exclude that height and CBP are also linked by biological pleiotropy. Further more advanced studies should be conducted to corroborate our findings. Prior studies demonstrating the vital role of SOX5 in normal vertebrate development are a reminder that measurements of human height used for GWAS may also reflect vertebral column development; if associations with CBP and height are connected via development of the vertebral column, it will be difficult to distinguish pleiotropy and causality using genetic studies alone.

22SNP-based heritability in the current study (8%) was considerably lower than estimates from twin studies (~40%). This is a common situation with modern methods of estimating heritability using genotype data, since such estimates reflect only one aspect of narrow-sense heritability captured by the additive genetic components of common variants, excluding the contributions of rare variants, non-additive effects, epistasis, or gene-environment interactions. Similar to what is seen with many other human traits, there was significant enrichment of SNP-based heritability of CBP for genetic regions that are conserved in mammals. Despite the modest heritability of CBP (and other self-reported traits), we found significant and large magnitude genetic correlations between CBP and other phenotypes that may be risk factors for CBP or consequences of CBP, such as depression-, osteoarthritis- and obesity-related traits (but not height). Future GWAS of CBP may benefit from taking these relationships into account, either as covariates, or in multivariate GWAS designs.

23A distinguishing feature of the current study as compared to many other GWAS is that the CBP phenotype examined represents a symptom, rather than a disease or a biomarker. Although successful GWAS of self-reported symptoms have been conducted which replicate associations seen with more specific disease phenotypes, our findings highlight potential challenges of GWAS of CBP: despite being one of the largest international studies of CBP ever conducted, our study detected only 3 significant associations with CBP. Still larger sample sizes will be needed in future discovery efforts using this phenotype, or different genetic approaches will be needed. A consequence of the nonspecific nature of the CBP phenotype is that, unlike other musculoskeletal phenotypes such as osteoarthritis, the tissue correlates optimal for conducting functional follow-up studies of findings from CBP GWAS are very unclear. Most animal models for back pain rely on specific mechanisms of inducing pain, such as injuries to the intervertebral disc, zygapophyseal (‘facet’) joint, dorsal root ganglion, or muscle. However, each of these mechanisms likely explain only a certain portion of back pain cases, and do not encompass the important psychosocial aspects of pain and pain reporting that are highly relevant in humans. Despite the importance of psychosocial factors, our meta-GWAS findings are a reminder that structural/anatomic factors involving spinal degeneration, such as disc herniation or osteoarthritis of spinal structures (e.g. facet joints), remain potentially important contributors to the CBP. While our study accounted for age by statistical adjustment, the meta-analysis design including multiple cohorts of older adults may have led to an overrepresentation of genetic variants associated with age-related conditions, such as osteoarthritis. Future GWAS of CBP may benefit from a broader age range of participants, stratification by back pain subtypes, simultaneously studying CBP and spinal degeneration/fracture phenotypes, and examination of interactions between genetic markers for spinal degeneration and markers for pain processing or axonal signaling (including DCC and netrin-1).

24Strengths of our study include its multicohort design and large sample size. A potential limitation of our study was heterogeneity of the CBP phenotype used, a consequence of pooling data from numerous cohorts using different definitions. Although this approach helped identify genetic associations shared across CBP subtypes, it might obscure associations pertinent to specific subtypes of back pain. As an example, we examined chronic back pain rather than chronic low back pain, since few of the included cohorts had available question items that isolated the low back region specifically. Given the high agreement between general back pain questions and low back-specific questions, and since mid/upper back pain without concurrent low back pain is uncommon, we expect that our results largely reflect genetic associations with low back pain. Despite phenotype heterogeneity, which would be expected to bias the study towards the null, we successfully identified several associations of statistical significance. Recent efforts to standardize CBP definitions may help to limit phenotype heterogeneity in future meta-GWAS of CBP. Another aspect of the phenotype used in our meta-analysis was that individuals with back pain of less than 3–6 months duration were included as controls. This was done deliberately so as to focus on back pain of chronic duration as the phenotype of interest. That said, GWAS examining back pain of any duration, or analyses excluding those with non-chronic back pain, may find different results. Another possible study limitation was lack of independence in the replication sample of UK Biobank participants from UKB2, given the same study base and methods between the UKB1 and UKB2 subcohorts. Our secondary analyses using linear mixed-effect models demonstrated similar SNP-CBP associations for our top hits when accounting for relatedness within UKB1, but the problem of relatedness across the two subcohorts (UKB1 vs. UKB2) remains. Finally, a limitation of this meta-analysis was that only autosomal variants were analysed, since some included cohorts did not analyze the X chromosome.

25In summary, this meta-analysis of GWAS of CBP identified novel genetic associations with CBP at SOX5, CCDC26/GSDMC, and DCC. Analysis of data from other GWAS and functional genomics experiments suggest possible pleiotropic effects of these loci on other traits including cartilage, osteoarthritis, lumbar disc degeneration, depression, and height/vertebral development, and possible causal effects on CBP mediated through height.

26Discovery 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. 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.

27There is no “gold-standard” definition for CBP. Consistent with the most commonly accepted clinical and research definitions for CBP, 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.

28Details 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, or the Haplotype Reference Consortium. 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 (S1 Table).

29We 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. 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).

30The 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−8 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).

31For 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, potential regulatory functions (RegulomeDB score), and effects on gene expression using expression quantitative trait loci (eQTLs) of different tissue types (GTExv6) (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. We used data from the Roadmap Epigenomics Project to evaluate whether the lead variants at each locus and those in LD (r2>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..

32Because 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, 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. 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.

33Given 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.

34Finally, we used LDsr of summary-level GWAS results from the discovery stage to estimate heritability due to common autosomal SNPs and genetic correlations. 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. We used stratified LDsr to partition heritability across functional categories of the genome, using methods described previously. 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.