Contribution of rare variation to degenerative orthopedic diseases

TL;DR

This study identifies rare genetic variations linked to degenerative orthopedic diseases using a large-scale analysis of genetic data.

Contribution

The study identifies new genes (TET2 and SMAD6) and confirms a known gene (LRP5) associated with rare variants in degenerative orthopedic diseases.

Findings

The LRP5 gene was confirmed as linked to osteoporosis with strong statistical significance.

TET2 was found to be associated with dorsalgia, and SMAD6 with spinal stenosis and other spondylopathies.

Only a small percentage of controls carried the identified rare qualifying variants.

Abstract

Degenerative orthopedic diseases (DODs) such as osteoporosis, osteoarthritis, spondylosis, spinal stenosis, and disc herniation are common disorders. Both common and rare genetic risk factors may contribute to DODs, but few large-scale whole-exome sequencing studies elucidating the contribution of rare variations to DODs have been published. The updated version of the Astra Zeneca portal (https://azphewas.com) was used to access gene collapsing analysis of rare variations for DODs.

One published UK Biobank portal: the updated Astra Zeneca portal based on whole genome sequencing (N = 484,111), was used to access gene collapsing analysis of rare qualifying variants (QVs) for fourteen DODs. A conservative threshold (p ≤ 5 × 10−8) was used to decrease the risk of spurious associations. Odds ratios (ORs) with 95 % confidence intervals (CIs) were estimated.

One previously osteoporosis-linked gene (LRP5) was identified (p-value = 3.0 × 10−13, OR = 1.55 (95%CI 1.38–1.73)). Two other significant genes were identified: the TET2 gene that was linked to dorsalgia (p-value = 8.3 × 10−9, OR = 1.54 (95%CI 1.33–1.77)) and the SMAD6 gene that was associated with spinal stenosis (p-value = 1.7 × 10−9, OR = 2.70 (95%CI 2.02–3.61)) and other spondylopathies (p-value = 9.0 × 10−11, OR = 2.41 (95%CI 1.89–3.07)). Among controls 1.10 % were carriers of LRP5 QVs, 0.29 % of controls carried TET2 QVs, and 0.20 %–0.21 % of controls carried SMAD6 QVs.

One established osteoporosis gene (LRP5) and two other DODs genes (TET2 and SMAD6) were identified using gene collapsing analyses. Only a small number of genes reached the chosen significance threshold under the current phenotype definitions and analytic framework.

Full text

Introduction

1

Degenerative orthopedic diseases (DOD) such as osteoarthritis, osteoporosis, and back pain due to spondylosis, spinal stenosis, and intervertebral disc degeneration are common and complex disorders where both genetic and acquired factors are involved [1,2]. Genome wide association studies (GWAS) have identified an increasing number of common variants associated with several DODs [[3], [4], [5]]. Exome sequencing may reveal causal rare variants for common complex disorders [[6], [7], [8], [9], [10], [11]]. However, most efforts for complex disorders have relied on GWAS, which focuses on common variants. Whole-exome sequencing (WES) or whole genome sequencing (WGS) studies can reveal the contribution of rare variants to common diseases such as DODs. However, only few WES studies of DODs have been published [[6], [7], [8]]. We used a published but updated UK biobank portal (https://azphewas.com/) [9] to access gene collapsing analysis of rare variation for DODs.

Methods

2

The Astra Zeneca portal published by Wang et al. originally reported the correlation between rare protein-coding variants and 17,361 binary phenotypes using WES data derived from 269,171 UK Biobank participants (https://azphewas.com/) [9]. The last version of the UK Biobank 500k WGS (v2) Public has been updated and is based on 484,111 whole-genome sequenced (WGS) individuals (https://azphewas.com/). The smaller Genebass portal (https://app.genebass.org/) that contains 394,841 UK biobank was not used [10]. All used data are publicly available, and ethical approval or informed consent was therefore not applicable. The gene-based collapsing approach, in which variants that satisfy specific criteria (qualifying variants, QV) are binned together as equivalent, is a simple but effective approach to identifying rare-variant contributions to disease [[9], [10], [11]]. The number of cases and controls with rare qualifying mutations in the gene are thereafter compared. Only the model with the lowest p-value of 12 tested models was selected from the Astra Zeneca portal (https://azphewas.com) [9]. In order not to obtain false positive associations, we used a more conservative p-value often used in GWAS studies (p ≤ 5 × 10^−8^) [[4], [5], [6], [7],[12], [13], [14]] instead of the commonly used p-value for WES studies (p-values <0.05/20,000 genes = 2.5 × 10^−6^) [[3], [4], [5],[9], [10], [11]]. In fact, the prespecified used threshold of 5 × 10^−8^ is close to the Bonferroni corrected p-value for WES studies (<0.05/20,000 genes/14 tested phenotypes/12 tested models = 1.5 × 10^−8^). The results are exploratory and should be interpreted under the chosen significance framework. We only show the union of the three-digit ICD-10 codes (International Classification of Diseases 10th Revision) [[9], [10], [11]]. Fourteen different DODs were included (Table 1): Osteoarthritis (M15-M19); back diseases (M47-M54), and osteoporosis with and without pathological fracture (M80, M81). Bioinformatic analysis was performed using Genecards (https://www.genecards.org) [12], OMIM (https://www.omim.org/) [13], and the GWAS catalog (https://www.ebi.ac.uk/gwas/) [14].Table 1. Results of gene collapsing analysis of rare variants for Degenerative orthopedic diseases (DODs) according to ICD-10 codes in the Astra Zeneca Portal (https://azphewas.com/) [9]. Union was used to define phenotypes for https://azphewas.com (updated version Dataset: UK Biobank 500k WGS (v2) Public). Only genes meeting the prespecified p ≤ 5 × 10^−8^ threshold are shown. Data was accessed 20251026. Odds ratios (ORs) are given with 95 % confidence intervals (CIs). For significant variants number of individuals and percentage with qualifying variants (QV) among cases and controls are given. Only the model with the lowest p-value of 12 tested models was selected from the Astra Zeneca portal (https://azphewas.com) [9].Table 1ICD-10 codesCases (QV n, %)Controls (QV n, %)Gene, p-value, OR (95%CI)M15 Polyarthrosis15,448122,880Not significantM16 Coxarthrosis27,444137,961Not significantM17 Gonarthrosis43,791155,123Not significantM18 Arthrosis of first carpometacarpal joint2848109,830Not significantM19 other arthrosis86,172132,789Not significantM47 spondylosis34,472138,266Not significantaM48 other spondylopathies16,513, (82, 0.50 %)148,889, (308, 0.21 %)SMAD6 9.0 × 10^−11^, OR = 2.41 (1.89–3.07)aM48.0 spinal stenosis9,939, (54, 0.54 %)154,052, (311, 0.20 %)SMAD6 1.7 × 10^−9^, OR = 2.70 (2.02–3.61)M50 Cervical disc disorders4781157,699Not significantM51 other intervertebral disc disorders26,819151,286Not significantM53 other dorsopathies not elsewhere classified6041119,393Not significantbM54 dorsalgia75,189, (333, 0.44 %)144,041, (416, 0.29 %)TET2 8.3 × 10^−9^, OR = 1.54 (1.33–1.77)M80 osteoporosis with pathological fracture2526113,109Not significantbM81 osteoporosis without pathological fracture24,886, (422, 1.70 %)97,889, (1,080, 1.10 %)LRP5 3.0 × 10^−13^, OR = 1.55 (1.38–1.73)Abbreviations: WGS = whole genome sequencing; ICD-10 = International Classification of Diseases 10th Revision; For Astra Zeneca portal (https://azphewas.com), the lowest p-value of the 12 tested models is shown [9].aUR: Model includes non-synonymous, predicted to be damaging (REVEL score ≥0.25), and Ultra-rare (minor allele frequency ≤0.00005 within the cohort and absent from gnomAD) (https://azphewas.com/) [9].bptvraredmg: Model includes non-synonymous, predicted to be damaging (REVEL score ≥0.25), PTVs (protein truncating variants) must be moderately rare (minor allele frequency ≤0.001 within the cohort and within gnomAD), and non-PTVs must be rare (minor allele frequency ≤0.00025 within the cohort and gnomAD global_raw ≤0.00005 and popmax ≤0.0005) (https://azphewas.com/) [9].

Results

3

In Table 1, only genes meeting the prespecified p ≤ 5 × 10^−8^ threshold are shown. Only the model with the lowest p-value of 12 tested models was selected from the Astra Zeneca portal (https://azphewas.com) [9]. One previously osteoporosis-linked genes (LRP5) were identified (Table 1) [[12], [13], [14]]. The only other two significant genes were the TET2 gene that was linked to dorsalgia (M54) and the SMAD6 gene that was associated with spinal stenosis (M48.0) and other spondylopathies (M48). (Table 1). No other DODs or genes were associated in the rare variant gene collapsing analysis (Table 1). Among controls 1.10 % were carriers of LRP5 QVs, 0.29 % of controls carried QVs for TET2, and 0.20 %–0.21 % of controls carried QVs for SMAD6 (Table 1).

No SMAD6 variants, i.e. rare or common, were linked to M48 or M48.0 in the variant analyses in the Astra Zeneca portal (https://azphewas.com) [9] or the GWAS catalog (https://www.ebi.ac.uk/gwas/) [14]. Nor were any rare or common TET2 variants linked to M54 in the variant analyses of the Astra Zeneca portal (https://azphewas.com) [9] or the GWAS catalog (https://www.ebi.ac.uk/gwas/) [14]. However, several LRP5 variants were linked to M81 in the variant analysis of the Astra Zeneca portal (https://azphewas.com) [9] and the GWAS catalog (https://www.ebi.ac.uk/gwas/) [14]. Thus, common variants are unlikely to contribute to the associations with SMAD6 and TET2 but may contribute to the association with the LRP gene.

The biological context and prior evidence as the result of bioinformatic search are presented in Table 2. The TET2 gene is linked to Myelodysplastic syndrome and Clonal hematopoiesis of indeterminate potential (CHIP) [12,13,15]. CHIP genes like TET2, DNMT3A, and ASXL1 have not previously been linked to DODs [[12], [13], [14]] Mutations in CHIP genes like TET2 may cause an increase in proinflammatory cytokines, resulting in chronic inflammation and immune imbalance and have been linked to many disorders [15]. The SMAD6 gene is linked to aortic valve disease, craniosynostosis, and radioulnar synostosis but not previously to DODs [[12], [13], [14]]. The SMAD6 gene encodes for a protein that belongs to the SMAD family of proteins, which are related to Drosophila ‘mothers against decapentaplegic' and C. elegans [12,13]. Among its related pathways are gene expression (transcription) and transforming growth factor beta (TGF-beta) receptor signaling in skeletal dysplasias [12,13]. SMAD6 is an inhibitory Smad (i-Smad) that negatively regulates signaling downstream of type I transforming growth factor-beta and acts as a mediator of TGF-beta and Bone Morphogenetic Protein anti-inflammatory activities [12,13].Table 2. Biological context and prior evidence. Bioinformatic search of the three genes associated with Degenerative orthopedic diseases (DODs) in the gene collapsing analysis of rare variants in the published UK Biobank portal (https://azphewas.com/ [9]. The GWAS catalog (https://www.ebi.ac.uk/gwas/) [14], OMIM (https://www.omim.org/) [13], and Genecards (https://www.genecards.org) [12] were searched for relevant diseases and pathways.Table 2. GenesGWAS catalog^12,14^Genecards/OMIM^12,13^Pathway analysis^12^SMAD6Bone tissue densityAortic valve diseaseGene expression (transcription)body heightCraniosynostosisTGF-beta receptor signaling in skeletal dysplasiasRadioulnar synostosisTET2body heightImmunodeficiency 75Myelodysplastic syndrome.MeiosisGene expression (transcription)LRP5Bone tissue densityOsteoporosis-pseudoglioma syndromeAlzheimer's disease and miRNA effects ncRNAs involved in Wnt signaling in hepatocellular carcinomabody heightOsteopetrosis, Autosomal dominantheel bone mineral densityBone mineral density variabilityEndosteal hyperostosisExudative vitreoretinopathyPolycystic liver disease with/without kidney cystsAbbreviations: GWAS = genome-wide association study; OMIM=Online Mendelian Inheritance in Man [13]. For a complete list of all GWAS associations, the reader is referred to the original bioinformatic resources [[12], [13], [14]].

Discussion

4

The updated version of the portal, UK Biobank 500k WGS (v2), based on 484,111 genomes [9] identified one well established osteoporosis gene (LRP5) and two other DODs genes (SMAD6 and TET2). However, the present study does not suggest that rare variation is a common contributor to common degenerative orthopedic diseases, though QVs in the SMAD6 gene are strong risk factors for spinal stenosis and other spondylopathies (Table 1). Thus, more common variants have been associated with DODs, especially osteoarthritis and osteoporosis [[3], [4], [5]]. There are no overt connections between the SMAD6 and TET2 genes and DODs but due to their involvement in multiple biological processes including regulation of inflammation, it is plausible that the UK biobank findings are correct [[12], [13], [14], [15]] (Table 2).

A strength of UK biobank is the large number of sequenced exomes, i.e. 484,111 in the last version (https://azphewas.com/) [9]. A limitation is the lack of validation studies in independent cohorts and more diverse populations concerning our findings, which will require large studies. It is, however, reassuring that one previously osteoporosis-linked gene (LRP5) [[12], [13], [14]] was identified. Other limitations are the lack of experimental evidence, and the use of a 3-digit ICD-10 code for the diagnosis of DODs in the UK Biobank without knowledge of the exact diagnostic procedures. Moreover, owing to the large number of comparisons, some associations might be spurious; we therefore applied a more conservative threshold (P ≤ 5 × 10^−8^) than the genome-wide threshold commonly used for WES studies (2.5 × 10^−6^) [[3], [4], [5],[9], [10], [11]]. On the other hand, potential candidates may therefore be discarded.

A limitation regarding TET2 is that no sensitivity analysis has been done with exclusion of established CHIP-driver mutations in TET2 to see if the association with dorsalgia remains. It is therefore possible that the finding is driven by age-related clonal hematopoiesis rather than a specific mechanism involving back pain.

The use of different models for the rare variant analysis in the portal (https://azphewas.com/) [9] is both a strength and a limitation. The use of 12 different models and selection of the model with the lowest p-value from the Astra Zeneca portal increases the chance of detecting more associations but also the risk for spurious associations. Another limitation is that rare variants can be highly population specific and findings derived from predominantly European datasets may not be generalizable to other ethnic groups. Some variants may also be limited to selected patients and families.

In conclusion, rare variations were significantly associated with DODs in the LRP5, TET2, and SMAD6 genes. Thus, rare variation contributes to DODs in the general population, though only a small number of genes reached the chosen significance threshold under the current phenotype definitions and analytic framework. Future rare variant analysis studies need to incorporate diverse populations or include selected patients and families with DODs.

Author contributions

All authors made substantial contributions to the conception and design of the work, and the interpretation of data. CAH and BZ drafted the original manuscript. All tother authors substantially contributed to the revision of the manuscript drafts. All authors have approved the submitted version of the manuscript and agreed to be accountable for any part of the work.

Funding

This work was supported by a grant awarded to Dr Bengt Zöller by ALF-funding from 10.13039/501100009780Region Skåne, and by the 10.13039/501100004359Swedish Research Council. The funders had no role in the study.

Conflicts of interest

None.