Sequence variation at 8q24.21 and risk of back pain

Introduction

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.

Methods

Results and discussion

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.

Figure 1. Pleiotropic effects of locus 8q24.21 on BP and related traits.

The BP-related traits and genes analyzed in our work are blue. The locus 8q24.21 is connected to associated traits (p-value < 5*10^-8) by green arrows. Possible relations between traits, as inferred from previous studies, are marked by grey arrows. Pairwise pleiotropic effects tested by SMR-HEIDI are denoted by a red cross if the test failed, and by a blue check if the test passed. A new detected pleiotropic effect on expression of FAM49B and BP based on SMR-HEIDI results is denoted by the blue arrow. LDH, lumbar disk herniation; SRDP, self-reported disk problems; HBMD, heel bone mineral density t-score.

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.

Conclusion

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.

Data availability