A multi-trait approach identified 7 novel genes for back pain.

The study design and the main results are presented in Figure 1 . Design of the study and the main results obtained in different steps. Venn diagram represents the distribution of phenotypes among individuals in the 500K UK Biobank data set. CBP, chronic back pain; GBA, gene-based association analysis; GWAS, genome-wide association study; IDD, intervertebral disk disorder; SGIT, shared genetic impact trait. The first step of the analysis was performed using a discovery sample from the UK Biobank. Initially, the GWAS summary statistics for the 3 BP-related phenotypes were obtained. They were used for the calculation of the multi-trait (shared genetic impact trait [SGIT]) summary statistics. Then we performed a gene-based association analysis (GBA) of the individual and multi-trait phenotypes. For all identified genes, conditional analysis (COJO) was applied. The second step of the analysis used a replication sample from the FinnGen database. We performed a replication study by applying the GBA to 32 genes identified in the first step and using GWAS summary statistics for 2 BP-related traits available on the FinnGen sample. A brief description of the methods is given below. More information can be found in the Supplementary Methods section, http://links.lww.com/PR9/A268 . Three BP-related conditions (CBP, dorsalgia, and IDD) were analyzed in individuals of White European ancestry from the UK Biobank 500k data set. We used genotyping and imputation data from the UK Biobank March 2018 data release. Genotypes of 19,397,568 SNPs, INFO >0.8, and MAF >5 × 10 −5 were analyzed. Genetic variants were annotated using the Ensembl Variant Effect Predictor (VEP). 18 We analyzed 3 variant annotations: protein coding (exons), protein noncoding (introns, 5′UTR, and 3′UTR), and nonsynonymous SNPs. The latter included transcript ablation, frameshift, stop gained, stop lost, start lost, transcript amplification, inframe insertion, inframe deletion, missense, and protein-altering variants. The numbers of the variants in each variant annotation are presented in Supplementary Table S1, http://links.lww.com/PR9/A268 . GWAS summary statistics were calculated using the fastGWA-GLMM tool, version 1.94.0 beta 10 and filters keeping variants with MAF ≥5 × 10 −5 and variant missingness rate ≤0.02. We included sex (data field 22001), age (data field 34), genotyping batch (data field 22828), and the first 10 genetic principal components provided by the UK Biobank as covariates in the association analysis to correct for fixed effects. To account for random effects mediated by relatedness, we used a sparse genomic relationship matrix (see Supplementary Methods for details, http://links.lww.com/PR9/A268 ). The summary statistics ( z -scores and effect sizes) were then used as input to 3 gene-based methods: SKAT-O, 12 PCA, 29 and ACAT-V. 16 These methods were implemented in the sumFREGAT R-package. 26 The results of these methods were combined by the aggregated Cauchy omnibus test, ACAT-O. 16 The significance level was defined as 2.5 × 10 −6 . For genes that passed significance threshold, gene-based analysis was repeated using conditional summary statistics obtained with the GCTA-COJO method. 30 , 31 Within the SHAHER framework, 27 a multi-trait phenotype, SGIT, was created by linearly combining the original traits. A GWAS of SGIT was then performed using an algorithm based on the summary statistics calculated for the original traits. SHAHER adjusts for sample overlap taking into account the phenotypic covariance between GWAS summary statistics. 27 The gene-based analysis was performed in the usual way. The SHAHER method was previously described in our studies 27 , 33 and is explained in detail in the Supplementary Methods, http://links.lww.com/PR9/A268 . GWAS summary statistics for phenotypes of “dorsalgia” and “other intervertebral disk disorders” of the FinnGen release v.9 11 were used to obtain gene-based P -values for statistically significant genes. Similar to the corresponding UK Biobank traits, the FinnGen phenotypes are electronic health record (EHR)-based and defined by ICD-10 codes (see https://r9.risteys.finngen.fi/endpoints/M13_DORSALGIA and https://r9.risteys.finngen.fi/endpoints/M13_INTERVERTEB for details). Chronic back pain summary statistics were unavailable because they were not included in the FinnGen release v.9. Summary statistics were filtered by INFO >0.8 and MAF >5 × 10 −5 . A replication threshold P -value was calculated as 0.05/(number of genes selected for replication). A gene was considered replicated if the P -value for either “dorsalgia” or “other intervertebral disk disorders” phenotype in the FinnGen sample was lower than the replication threshold.

We observed 86 genome-wide significant gene-based signals across all trait-annotation combinations, with 54 unique genes detected in total: 51 for CBP, 3 for IDD, and 2 for dorsalgia (2 genes associated with both CBP and IDD, see Supplementary Table S2, http://links.lww.com/PR9/A268 ). Full results of the gene-based analysis are available in the ZENODO database ( https://doi.org/10.5281/zenodo.8118630 ). QQ-plots for gene-based analyses are shown in Supplementary Figure S1, http://links.lww.com/PR9/A268 . Manhattan and QQ plots for single-point association analyses are shown in Supplementary Figure S2, http://links.lww.com/PR9/A268 . Using the phenotypic and genetic correlations between the original traits and their SNP-based heritability (Supplementary Figure S3, http://links.lww.com/PR9/A268 , first 3 columns and rows), we estimated the coefficients of optimal linear combination of the original traits to build SGIT as 0.69, 0.44, and 0.31 for CBP, dorsalgia, and IDD, respectively. Then, we calculated the phenotypic and genetic correlations between the original traits and SGIT and estimated the SNP-based heritability of SGIT (Supplementary Figure S3, http://links.lww.com/PR9/A268 , fourth column and row). The phenotypic correlations between SGIT and the original traits are stronger than those between the original traits. All genetic correlations between SGIT and the original traits are no less than 0.93. The SNP-based heritability of SGIT is higher than that of the individual traits. We observed 105 signals in 65 genes (Supplementary Table S2, http://links.lww.com/PR9/A268 ). We tested a total of 191 signals in 87 genes identified using the original traits and SGIT. After conditional analysis, 15 of 80 signals obtained for CBP retained significance. For IDD, 1 of 4 signals passed conditional analysis. For SGIT, 25 of 105 signals retained significance after conditional analysis. For dorsalgia, neither of the 2 genome-wide significant signals passed COJO (Supplementary Table S2, http://links.lww.com/PR9/A268 ). In total, 32 genes were selected for replication. We tested the replication of 32 genes left after the conditional analysis using FinnGen summary statistics obtained for IDD and dorsalgia. The threshold P -value level for replication was calculated as 0.05/32 = 0.00156. A total of 16 genes were replicated (Table 1 ). All genes were replicated using noncoding variant annotation. The CHST3 gene identified for IDD passed the replication threshold. Among 15 genes significantly associated with CBP after conditional analysis, 2 genes, C8orf34 and RHOA , were replicated. Among 24 genes identified for SGIT, 13 genes were replicated. Six of them, DCC , ILRUN , MAML3 , MIPOL1 , SOX5 , and SPOCK2 , were also replicated for CBP (Table 1 ). Manhattan plots for gene-based analysis using noncoding, coding, and nonsynonymous variants are shown in Figure 2 , Supplementary Figures S4 and S5, http://links.lww.com/PR9/A268 , respectively. Replication results on FinnGen data. Genes that passed the replication threshold are shown in bold. “cod”: protein coding SNPs; “ncod”: protein noncoding SNPs; “nsyn”: nonsynonymous SNPs. CBP, chronic back pain; IDD, intervertebral disk disease; SGIT, shared genetic impact trait. Manhattan plots for gene-based association analysis using noncoding variants. Blue and red triangles indicate −log10 ( P value) for analyzed genes before and after conditional analysis, respectively. Replicated genes are labeled. Dotted line indicates the significance level (2.5 × 10 −6 ). For 16 replicated genes, we compared the P -values obtained on SGIT and original BP-related phenotypes (Table 2 ). An overlap between association analysis results of SGIT and original traits is shown in Figure 3 . Results of association analysis of original traits for genes identified on shared genetic impact trait. CBP, chronic back pain; IDD, intervertebral disk disease; SGIT, shared genetic impact trait. Overlap of association analysis results for multi-trait and original traits. CBP, chronic back pain; IDD, intervertebral disk disorder; SGIT, shared genetic impact trait.

We have performed gene-based analysis of 3 BP-related phenotypes and their multi-trait combination using the UK Biobank data. We identified and replicated 16 genes, 13 of which were detected using multi-trait phenotype SGIT. Seven genes identified using SGIT were undetected using the original BP-related phenotypes. In recent years, several GWAS of BP-related phenotypes have been performed using the UK Biobank data. 6 , 24 However, only a few independent loci were identified and replicated. 14 The success of our study in detecting BP-related genes can be explained by the new methodology used for analysis. In our study, 2 approaches have been newly applied for BP-related phenotypes: gene-based and multi-trait association analyses. The gene-based association analysis of BP has been previously performed only for the analysis of rare variants from exome sequenced data 1 , 33 but not imputed genotypes. In addition, a multi-trait approach was used to solve the following problem. According to a contemporary model, there is a complex interplay between biological, psychological, and social factors that contribute to back pain. 7 The identification of novel genetic loci is severely hampered by the absence of a consensus definition of BP. To overcome these obstacles and improve our comprehension of the genetic impact on BP, we used a novel approach. Three distinct phenotypes related to back pain were taken into consideration: 2 were based on EHR-based codes and 1 was self-reported. Dorsalgia and IDD are 2 EHR-based back pain–related phenotypes that are different from self-reported CBP in that they are based on longitudinal health care records, frequently spanning many years, rather than point prevalence from a single survey. They reflect back pain that is severe enough to warrant seeking medical attention. Nevertheless, because of the roughly 90% genetic correlations among these phenotypes, the genetic nature of all of them is common. However, it does not mean that all 3 BP-related phenotypes will manifest in the majority of the patients. The discordance between phenotypes is explained by the relatively low heritability (approximately 50%) of each phenotype. Because of this, only some patients manifested more than 1 BP-related phenotype (Fig. 1 ). We explored shared pathways among different pain phenotypes using multi-trait methodologies, which can help reduce heterogeneity. 27 , 28 We used SHAHER methodology 27 for multi-trait analysis. SHAHER constructs a multi-trait phenotype as a linear combination of original traits, where the coefficients are estimated to maximize the heritability of the multi-trait phenotype. As a result, we identified 7 genes associated with SGIT, which did not show a significant association with the original traits. Additional 6 genes were associated with both SGIT and CBP. For the majority of genes associated with SGIT, the P -values obtained on the original BP-related traits were rather low (Table 2 ). These results confirm the effectiveness of the chosen strategy. Earlier, we applied both methodologies to analyze the associations between BP and rare genetic variants using exome sequenced data. 33 We identified a new FSCN3 gene associated with BP-related multi-trait phenotype because of its loss of function (LOF) variants. These variants changing protein structure can be identified when exome sequenced data is used. Imputed genotypes have better coverage of noncoding variants. All genes identified and replicated in the current study showed an association because of noncoding variants (Table 1 ). Therefore, association analysis using imputed and sequenced genotypes provided different knowledge about the genetic architecture of BP. The highest association signals were detected for 2 genes, SOX5 ( P = 1.8 × 10 −11 for SGIT and P = 7.5 × 10 −13 for CBP) and SPOCK2 ( P = 8.1 × 10 −13 for SGIT and P = 2.7 × 10 −10 for CBP). The SOX5 gene codes transcription factors that play essential roles in chondrocyte differentiation 32 ; this gene has been previously detected in Reference 3 and by our group as associated with BP-related phenotypes. 3 , 6 , 24 The association of SPOCK2 was described in Reference 6 . It is interesting that SPOCK2 is located close to the CHST3 gene, which codes for an enzyme that catalyzes proteoglycan sulfating and has been identified as a susceptibility gene for lumbar disk degeneration. 23 The association of CHST3 with IDD and dorsalgia was described originally in a study of Asian participants. 3 Moreover, this gene is considered as most probably causal for back pain in this region (cumulative evidence L2G scores for CHST3 and SPOCK2 are 0.81 and 0.23, respectively, https://genetics.opentargets.org/study-locus/NEALE2_6159_4/10_72001257_A_G ). 19 We believe that a significant result for the SPOCK2 gene indicates an inability of association analysis to resolve the causal gene in the region, especially in cases with a strong LD block. In such instances, using GWAS results from other ancestries has the potential to provide assistance. For the additional 6 genes identified: BSN , DCC , ILRUN , SMAD3 , COL11A , and C8orf34 , the association with BP-related phenotypes has already been shown, 3 , 24 and our results can be considered as verification of the impact of these genes on BP-related traits. For the other 7 genes, no direct association with back pain phenotypes has been shown yet. However, 3 of them were located in known loci, but have not been prioritized: MIPOL1/SLC25A21 , PTPRC / MIR181A1HG , and RHOA/BSN . 3 Current data do not provide enough information to confidently choose a single candidate gene within these loci. For PTPRC , the association with inflammatory bowel disease including Crohn disease has been shown. 15 The inflammatory conditions have significant genetic correlations with back and multisite chronic pain. 22 MAML3 has been shown to contribute to several pathways significantly associated with CBP. 4 This gene enables transcription coactivator activity; it is involved in the Notch signaling pathway and positive regulation of transcription by RNA polymerase II. For 3 further genes: RERG , JADE2 , and MLLT10 , association has been found in a multi-trait analysis of endometriosis and multisite chronic pain. 22 In our study, we observed these genes to be associated with SGIT but not the other phenotypes. As they seem to affect many comorbid pain conditions, this is in line with the multi-trait nature of SGIT. This conclusion is supported by the fact that some genes associated with SGIT are involved in chondrogenesis and inflammatory processes, which are common for BP conditions. In accordance with the GWAS catalog, many genes are associated with the neurological processes (cognitive ability, educational attainment, attention deficit, intelligence, depression, and others), immunological and inflammatory processes (Crohn disease, ulcerative colitis, asthma, osteoarthritis, ankylosing spondylitis), and bone mineral density. All these processes can be involved in BP. 3 , 20 , 24 However, the exact role of many identified genes in BP is unknown and needs to be clarified. Our study has several limitations. First, we used only White Europeans because the linear combination coefficients for SGIT cannot be estimated using GWAS of small sample sizes. Therefore, not all our findings and conclusions could be applied to other ethnic groups. Second, we replicated the association signals using the freely available FinnGen summary statistics for dorsalgia and IDD. The summary statistics for CBP were absent in this data set and were not included in replication. We replicated 16 of the 32 genes that passed conditional analysis. It is possible that some of the nonreplicated genes will be replicated later, when summary statistics obtained from other large projects will be accessible for replication.

We conducted a wide gene-based association analysis of BP-related phenotypes using imputed genotypes from a UK Biobank project and applying the technique of multi-trait analysis. We identified and replicated 16 genes significantly associated with BP-related traits. Thirteen genes were detected on multi-trait phenotype that is in accordance with high genetic correlation between BP-related traits. Some genes have been previously described as associated with BP or with the genetically correlated traits or as included in pathways associated with BP.

Back pain (BP) is an important contributor to global disability. 9 , 17 Chronic back pain (CBP), defined as self-reported back pain lasting more than 3 months, is the most debilitating form of back pain. It is prevalent in approximately 10%, depending on the population and the definition used. 17 Chronic back pain is a complex heritable trait with heritability ranging from 40% to 68%. 8 Genes associated with CBP have been identified in several genome-wide association studies (GWAS) using data from national biobanks. 3 , 6 , 24 , 25 , 28 Recent advances in BP genetics research have been made by meta-analysis of other BP-related phenotypes such as dorsalgia and intervertebral disk disease (IDD). 24 A total of 41 variants at 33 loci were identified by the authors. However, the liability scale heritability of BP estimated from the SNPs is only 13%, which is significantly lower than the 40% heritability of BP expected from classic twin studies. 21 Several explanations have been proposed to explain this lack of heritability. One of them is the need for more powerful methods of association analysis. 5 In this study, we used a novel multi-trait approach to merge 3 BP-related phenotypes: CBP, dorsalgia, and IDD, which were previously used in BP studies. Dorsalgia is a phenotype based on electronic health records (EHRs) and is defined by clinical diagnostic codes for back and neck pain, reflecting mainly the former. Intervertebral disk disease is another phenotype related to BP and defined by EHR-based codes for intervertebral disk degeneration. These codes are typically used when BP is present, but not necessarily. 2 , 3 There are high genetic correlations between BP-related phenotypes, estimated to be in the range of 87% to 92%, 3 , 33 allowing the use of methods that maximize multi-trait heritability. The aim of this approach was to identify a common genetic background for a number of traits that are genetically correlated. We also used gene-based association analysis, which is a robust method because it estimates the combined effect of a set of variants within a gene rather than analyzing each variant individually. 13 This method enables the identification of genes whose protein-coding mutations alter the structure of the corresponding proteins or whose intragenic mutations affect gene expression. Recently, we applied multi-trait and gene-based approaches to analyze the association between BP and rare exome-sequenced variants and identified a new FSCN3 gene. 3 , 33 In this study, our aim was to continue the analysis of the genes responsible for differences in the risk of BP between individuals using multi-trait gene-based association analysis. Previously, we did this by testing the effect of rare and ultra-rare variants of sequenced genotypes in the UK Biobank 200k dataset. 33 Here we analyze common variants using imputed genotypes for the larger UK Biobank dataset containing 500k individuals. Imputed variants are predominantly located in the regulatory intron regions of genes, whereas exome-sequenced variants are predominantly protein coding.

Supplemental digital content associated with this article can be found online at http://links.lww.com/PR9/A268 .

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.