Back pain (BP) is a common debilitating condition with a lifetime prevalence of 40% and a major socioeconomic impact ¹ . According to the Global Burden of Disease 2016 study, it leads the list of disabling conditions in many parts of the world ² . The greatest biological risk for episodes of severe BP is thought to be lumbar disc degeneration (LDD) ³ . Rather than a discrete disease, LDD is now considered a continuous lifelong aging process that usually starts in the third decade of life but may also be seen in children ⁴ . It is highly heritable ⁵ , as well as being influenced by smoking ⁶ , age ⁷ and body mass index ⁸ . The overall burden of lumbar disc degeneration on magnetic resonance imaging is a strong predictor of BP episodes ^{9– 11} . There is a clear genetic predisposition to both BP and LDD with estimates of heritability in the range of 30%–70% ^{12– 15} . Also, these traits have been shown to share genetic risk factors ¹⁶ . However, little is known about the genetic architecture of BP and only three associated genetic loci have been identified ^{17, 18} .

This work was focused on the locus at chromosome 8 (8q24.21), which was identified in our previous BP studies ^{17, 18} . The same locus has been reported associated with having lumbar disc herniation (LDH) in Icelanders undergoing surgery for sciatica ¹⁹ . Bjornsdottir et al. additionally reported an association of this region to height, but concluded that the region does not influence LDH through height. They noted that the most highly associated SNP, rs6651255, was linked to expression of the nearby gasdermin ( GSDMC) gene in GTEx data ²⁰ . Further, the immune-related function of gasdermins implied that the inflammatory and immunological consequences of disk herniation may be linked to sciatica symptom severity.

However, closely located SNPs associated with different traits do not always reflect pleiotropy ²¹ . Furthermore, in the case of eQTL studies, causal hypotheses generated through single-variant look up are highly likely to be false positives ²¹ . Consequently, a variety of co-localization methods have been introduced ^{22– 25} to distinguish the two common scenarios: pleiotropy (the same functional genetic variant is associated with both traits) and linkage (two different functional genetic variants, each associated with a different trait, that are in linkage disequilibrium).

The aim of this study was to (1) verify whether the same functional variant at 8q24.21 is responsible for BP and LDH and (2) to study the effects of this variation on related phenotypes. As a result, we offer an alternative explanation to the assertion by Bjornsdottir et al. ¹⁷ regarding the role of immunity and inflammation in generating symptoms of sciatica and BP due to variation in expression at the locus in between genes CCDC26 and GSDMC.

We had previously performed the largest-to-date GWAS for chronic BP using the UK Biobank and CHARGE consortium cohorts, with total of 119384 cases and 334478 controls ¹⁷ . Cases were defined as those having BP of duration >3-6 months; controls included those without back pain or with back pain of <3-6 months. We found variants at 8q41.21 locus that showed a genome-wide significant association (OR=1.05, p = 4.4 × 10 ^-13, after genomic control correction) ¹⁷ . The location chr8:130,717,716 is tagged by intergenic variant rs10956487 which is in perfect LD with rs6651255 (r ²=1 according to 1000 Genomes EUR samples), the lead SNP in the Icelandic study of LDH by Bjornsdottir et al. ¹⁹ . The strongest association signal for BP is situated near GSDMC gene; yet another protein-coding gene, FAM49B, is located within a 500 kbp window. The summary of the previously reported results and our findings for locus 8q24.21 are shown in Figure 1.

We used the SMR-HEIDI approach ²⁴ to test for potential pleiotropic effects of the 8q24.21 region (rs6651255) on BP and other complex traits including gene expression levels in relevant tissues. SMR analysis provides evidence for pleiotropy, but is unable to determine whether both traits are influenced by the same underlying causal polymorphism. The latter question was addressed by the HEIDI test which distinguishes pleiotropy from linkage. For SMR-HEIDI testing we used as the target trait the BP GWAS previously performed in all genetically confirmed white British UK Biobank participants ¹⁷ .

First of all, we checked relationships between BP and LDH. We found strong evidence that co-association of these traits to 8q24.21 is due to the pleiotropic action of the same functional variant (P _SMR=2.8 × 10 ^-7 and P _HEIDI = 0.87; extended data ²⁹ ). We further applied SMR-HEIDI to test for pleiotropic effects of the region in the collection of 2,453 traits available from the UK biobank study and GeneAtlas. We found that the BP risk allele of rs6651255 was negatively associated with heel bone mineral density t-score and positively associated with self-reported prolapsed or slipped disc (both FDR _SMR<0.01, underlying data ³⁰ and extended data ²⁹ ) and also strongly associated with standing height (underlying data ³⁰ and extended data ²⁹ ) and several measures of fat-free mass (all FDR _SMR<0.001). However, HEIDI tests showed strong evidence against pleiotropic effects of the same functional variant in this region for BP and these anthropometric traits (all P _HEIDI < 0.001), while associations with heel bone mineral density t-score and self-reported prolapsed or slipped disc were likely the result of pleiotropy, controlled by the same functional variant (all P _HEIDI > 0.05). This suggests that genetic variants influencing structural spine features – the vertebral body and intervertebral disc – provide the link to BP symptoms, supporting previous work in twins highlighting shared genetic influence between spine degeneration and BP ^{16, 31} .

The unanticipated results were obtained in the study of expression data. We have previously ¹⁷ demonstrated that variants in LD with the lead SNP at 8q24.1 were located in potential regulatory regions in chondrocytes and osteoblasts (Roadmap Epigenomics Consortium ³² ). In the current study, we tested expression of genes in the 8q24.21 region in all available tissues considered potentially relevant to development of pain (underlying data ³³ ), including musculoskeletal, brain and nerve tissue (from GTEx, ²⁰ ). Given the suggestion by Bjornsdottir et al. that the locus may act through regulation of immune function ¹⁹ , we also investigated expression in whole blood using the blood eQTL database ²⁶ and five immune cell subsets from CEDAR ³⁴ . Similar to the Bjornsdottir et al. study ¹⁹ , we observed that rs6651255 is co-associated with expression of GSDMC in skeletal muscle (FDR _SMR=0.01, P _SMR=0.0001, underlying data ³⁰ ), with the BP-associated C allele of rs6651255 being associated with increased expression of the gene. However, HEIDI results show this overlap to be unlikely due to pleiotropic action of the same functional variant (P _HEIDI = 3×10 ^-7, underlying data ³⁰ and extended data ²⁹ ). Our findings suggest that GSDMC is not the gene through which the functional variant acts to influence BP. The next strongest co-association was with increased expression of FAM49B in brain anterior cingulate cortex (P _SMR=0.004, underlying data ³⁰ and extended data ²⁹ ) and decreased expression in cytotoxic T (CD8+) cells (P _SMR=0.02, underlying data ³⁰ and extended data ²⁹ ), although these associations were only nominally significant (p<0.05), yet not statistically significant after multiple testing correction (FDR _SMR=0.18 and 0.28, respectively). All corresponding P _HEIDI were > 0.05, indicating these co-associations are likely due to pleiotropy. The product of FAM49B was found to be expressed in multiple tissues including bone, lymph node and spleen; according to Proteomics DB ³⁵ , it has highest expression in lymphocytes, and was implicated as a key negative regulator of actin dynamics and T-cell activation ³⁶ . The precise mechanism of which gene is involved in the development of BP is difficult to infer due to lack of eQTL data available for relevant tissues, such as intervertebral disc.

We hypothesize that the locus 8q24.21 exhibits its influence on BP phenotypes via disk and bone remodeling, and not through influence on anthropometrics or pain processing. Our findings suggest that FAM49B and not GSDMC is the most likely mechanism through which the functional variant leads to BP.