Osteoarthritis is one of the most rapidly increasing health conditions globally, and among the leading causes of disability and pain1. The global burden of osteoarthritis has reached a staggering 595 million individuals, representing a notable 132% increase in prevalence since 19902. The total number of patients with osteoarthritis has been estimated to reach 1 billion worldwide by 20502. Despite the enormous societal and public health burden of osteoarthritis, no effective disease-modifying treatments exist. It is therefore imperative to enhance our understanding of the biological processes leading to disease development to accelerate translation.
Osteoarthritis is a complex disease, caused by an interplay between environmental and genetic risk factors. Previous genome-wide association studies (GWASs) have led to the identification of around 150 risk variants, mediated through effector genes involved in various pathways3. Here we conducted a large-scale GWAS meta-analysis across 1,962,069 individuals, achieving a 2.64-fold increase in effective sample size compared with the next largest GWAS3. We combine the genetic findings with functional genomics evidence from osteoarthritis-relevant tissues and identify effector genes that converge on key biological processes underpinning disease development, generating insights into targets for focused therapeutic interventions.
We have performed a large multi-ancestry GWAS meta-analysis for osteoarthritis, combining 87 datasets across 489,975 cases and 1,472,094 controls, with an effective sample size of 1,470,467 individuals (Methods and Supplementary Table 1). It includes 87.31% individuals of European (EUR) ancestry, 7.09% East Asian (EAS) ancestry, 3.08% African American (AFR) ancestry, 1.09% South Asian (SAS) ancestry, 0.91% Hispanic (HIS) ancestry and 0.53% with mixed ancestry (ADM) (Supplementary Tables 1 and 2). In addition to osteoarthritis at any joint as an overarching disease phenotype, we performed joint-specific GWAS meta-analyses on the basis of the joint affected (Methods).
We identified 962 independent osteoarthritis associations at the study-wide significance threshold of P ≤ 1.3 × 10−8 (175 for osteoarthritis at any site, 151 for hip osteoarthritis, 146 for knee osteoarthritis, 131 for hip and/or knee osteoarthritis, 4 for spine osteoarthritis, 14 for hand osteoarthritis, 7 for finger osteoarthritis, 5 for thumb osteoarthritis, 136 for total hip replacement, 92 for total knee replacement and 101 for total joint replacement) (Fig. 1, Supplementary Figs. 1–3 and Supplementary Table 3), some of which overlap across phenotypes. The majority of these (513 out of 962) are conditionally independent of any previously reported risk variant for any osteoarthritis phenotype (Supplementary Tables 3 and 4). Of the 962 variants, 339 are unique and conditionally independent across all osteoarthritis phenotypes (236 newly reported here) (Methods).
Meta-analysis-based odds ratios of the 962 index variants as a function of their risk-allele frequency, and phenotypic variance explained (VEP) for each variant indicated by the size of each circle. Each colour corresponds to an osteoarthritis phenotype: osteoarthritis at any site (ALLOA), hip osteoarthritis (HIP), knee osteoarthritis (KNEE), hip and/or knee osteoarthritis (HIPKNEE), spine osteoarthritis (SPINE), hand osteoarthritis (HAND), finger osteoarthritis (FINGER), thumb osteoarthritis (THUMB), total hip replacement (THR), total knee replacement (TKR) and total hip and/or knee replacement (total joint replacement, TJR).
The 962 independently associated variants map to 286 genomic loci (176 newly reported here). Of the 110 previously reported loci, 44 have a newly reported, independent osteoarthritis-associated variant (Methods and Supplementary Tables 3 and 4). Most loci (86%) contain a single independent signal, with the remainder encompassing between 2 and 5 independent signals per locus. 95% of the associated variants have a minor allele frequency (MAF) of ≥5% with small to modest effects (odds ratios (OR), 1.016–1.186). Forty-nine signals are driven by low-frequency variants (MAF, 1–5%; OR, 1.044–1.279) (Fig. 1 and Supplementary Fig. 1). We performed GWAS meta-analysis within four ancestry groups (EAS, AFR, SAS, HIS) and did not detect ancestry-specific study-wide significant associations. We also did not find any additional signals when restricting the GWAS meta-analysis to studies in which osteoarthritis had been defined based on imaging. We find high correlation between associations when comparing GWAS with and without the inclusion of self-reported osteoarthritis (Methods, Supplementary Figs. 4–6 and Supplementary Note).
In addition to the 339 unique signals from the main analyses, we find 3 newly reported female-specific associations and 1 male-specific association with significant differences in effect size between sexes (Phet < 0.0125) that did not reach genome-wide significance in the combined sex analysis (Methods and Supplementary Tables 5 and 6).
We evaluated the predictive potential of genetic risk scores (GRSs) in independent datasets (Methods). For the osteoarthritis phenotypes tested, no analysis reached an area under the receiver operating characteristic curve (AUC) over 80%. The best-performing GRS was obtained for hip osteoarthritis (AUC, 58.6%) (Supplementary Table 7). We found that including body–mass index (BMI) in the GRS model led to improvements in prediction (for example, hip osteoarthritis including BMI AUC, 66%).
To determine whether early development of skeletal tissues contributes to the risk of osteoarthritis later in life, we investigated the enrichment of GWAS signals in cell types associated with skeletal development through functional GWAS (fGWAS). We performed the analysis for 30 different cell types using single-cell multiomics data (ATAC and RNA-seq) from the human skeletal development atlas4, spanning 5–11 weeks after conception (Fig. 2, Methods and Supplementary Table 8).
fGWAS enrichment for osteoarthritis in 30 cell states of the skeletal development atlas. Significance (FDR < 0.1) and effect size (log-transformed OR, log[OR]) are indicated by colour and dot size, respectively. InterzoneChon, interzone chondrocytes; PAX7high Chon, PAX7-expressing chondrocytes; ChondroPro1, chondrocyte progenitors; CyclingChon, chondrocytes with high cell cycle activity; ArticularChon1, articular chondrocytes with high TRPV4 and VEGFA expression; ArticularChon2, articular chondrocytes with high EPYC and low SOX9 expression; DLK1high Chon, DLK1-expressing chondrocytes; HypertrophicChon, hypertrophic chondrocytes; MaturingChon, maturing chondrocytes; LimbMes, early limb mesenchyme cells; Perichondrium, perichondrial osteoblast progenitors; MatureOsteocyte, osteocytes; FibroPRO1/2, fibroblast progenitors; SynFIB, synovial fibroblasts; DermFIB1/2, dermal fibroblasts; TENO, tenocytes; PAX7 Myo, PAX7-expressing myocytes; MYH3 Myo, MYH3 expressing myocytes; PERI, pericytes; PerineuralFIB, perineural fibroblasts; HIC1 Mes, HIC1-expressing mesenchymal cells.+++
In the chondrogenesis lineage, we observed significant enrichment (false discovery rate (FDR) < 0.1) for mature, hypertrophic, articular and DLK1-expressing chondrocytes for all of the tested osteoarthritis phenotypes, consistent with cartilage being the primary affected tissue. Chondrocytes with high cell cycle activity were also enriched for all phenotypes except for finger osteoarthritis. Moreover, more immature cell types including chondrocyte progenitors and early GDF5 expressing interzone chondrocytes were enriched for total hip replacement, and chondrocyte progenitors were also enriched for hip osteoarthritis. In the osteogenesis lineages, we observed significant enrichment for mature osteocytes (total hip replacement, hip osteoarthritis and finger osteoarthritis), osteoblast (total hip replacement and hip osteoarthritis) and perichondrium (hip osteoarthritis).
The osteoblast enrichment associated with hip and finger osteoarthritis may be linked to bone morphology, as structural abnormalities in femoral head formation can lead to irregular joint surfaces or improper joint congruity, increasing the risk of mechanical overloading, and contributing to osteoarthritis development. Geometric parameters of the hip are known to be associated with osteoarthritis5,6, and developmental dysplasia of the hip often leads to osteoarthritis, with research showing shared genetic risk factors between the two conditions, including associations with GDF5 and COL11A17,8. Finger-length patterns in combination with elevated androgen levels during development have also been linked with osteoarthritis9. The fGWAS results therefore suggest a role of bone development in the pathogenesis of hip and finger osteoarthritis manifesting in later stages of life and implicates particular transcriptomic and epigenetic cell states.
We find enrichment in total hip replacement and hip osteoarthritis genetic associations with tenocytes. Tendons are vital to the transmission of force and stabilization of the musculoskeletal system. Hip tendon samples from patients with osteoarthritis demonstrate a greater degree of fibrosis, non-collagenous change and calcium deposition in the extracellular matrix (ECM) compared with samples from patients with femoral neck fractures10, consistent with periarticular tendinopathy. Similar tendinopathy is found at other osteoarthritis-susceptible joints11,12. Our findings indicate that tendon development is also associated with hip osteoarthritis and is more likely related to late-stage osteoarthritis, suggesting that the developmental biology of secondary stabilizers of the joint contributes to the causal pathway in osteoarthritis.
To identify potential causal variants at the associated loci, we created, at each signal, a set of variants that are predicted with 95% probability to include a causal variant, called credible sets (Methods and Supplementary Note). The number of variants in a credible set ranged from 1 to 247 (mean 23 variants) with 75 credible sets containing a single variant and 149 credible sets containing less than 3 variants (Supplementary Tables 9 and 10). A total of 328 credible sets mapped entirely within the transcript of a single gene, strongly indicating that gene as causal. Most credible-set variants were predicted to be non-coding (57% were intronic and 17% intergenic). In total, 81 coding credible-set variants were missense, 1 was a stop gain variant (in VIT) and 1 was a splice acceptor variant. On the basis of 3D chromatin interaction data that we generated in primary osteoarthritis chondrocytes (Methods), 187 credible-set variants overlap promoters, 2,149 overlap enhancers and 814 reside within an enhancer that loops to a promoter. We performed transcription factor enrichment analysis (Methods) and identified 1,585 credible-set variants that both reside within gene regulatory regions and affect a transcription-factor-binding motif in osteoblast or chondrogenic cells (344 unique transcription factors; Supplementary Tables 11 and 12, Supplementary Fig. 7 and Supplementary Note).
To identify genes that are very likely to be causal for osteoarthritis (effector genes), we integrated data across 24 orthogonal lines of evidence to score each of the 8,785 genes residing within the 286 genomic risk loci (Methods, Extended Data Fig. 1 and Supplementary Tables 13–19). We identified 700 unique effector genes with a score of ≥3, mapping to over 88% of loci (Supplementary Table 20). We find that 70 loci contain a single effector gene, while the majority (70%) contain more than one gene with at least three orthogonal lines of evidence pointing to its involvement. The highest-scoring effector gene, with 11 lines of evidence in support of its involvement, is ALDH1A2, a gene previously implicated in osteoarthritis13.
We found that mouse and human musculoskeletal and pain phenotypes, chondrocyte HiC and differential chondrocyte methylation are the lines of evidence with relatively higher information contributions (Methods, Supplementary Tables 13, 21 and 22, Supplementary Fig. 8, Extended Data Fig. 2 and Supplementary Note).
We assessed the association between loss of function (LOF) variants in the 700 effector genes and osteoarthritis using gene burden tests. To this end, we aggregated the association of all rare LOF variants in these genes (<2% frequency total) (Methods and Supplementary Tables 23–25) and identified nine study-wide significant associations (P < 7.1 × 10−5) with 5 genes (ADAMTSL3, VIT, COL27A1, IL11 and PMVK), of which the burdens of ADAMTSL3 and VIT on hip osteoarthritis and total hip replacement are genome-wide significant (P < 2.5 × 10−6). The risk of disease was increased for LOF variants in these genes. When we incorporated missense (MIS) in addition to LOF variants (LOF + MIS) in the burden tests, we identified ADAMTSL3, VIT, IL11, THBS3, ADAMTS6, SPRY2 and COLGALT2 associated with osteoarthritis, of which association of ADAMTSL3 with hip osteoarthritis and IL11 with total hip replacement are genome-wide significant. LOF + MIS variants in ADAMTS6, SPRY2 and COLGALT2 are protective against osteoarthritis, whereas aggregation of these variants in ADAMTSL3, VIT, IL11 and THBS3 confer risk of osteoarthritis. The direction of effects was consistent in both models for all effector genes. Common non-coding sequence variants associated with osteoarthritis phenotypes present concordant directions of effect with gene-burden association results of genes in their vicinity, with the exception of variants near THBS3 and PMVK; these two genes are at the same locus (around 300 kb apart). Notably, none of the above burden associations are driven by a single variant in any of the cohorts (Supplementary Table 25).
We found LOF burdens for genes at the same loci as those identified in the common variant analysis for the same phenotypes and for different phenotypes (for example, ADAMTSL3 and total hip replacement, PMVK and knee osteoarthritis, and SPRY2 and hand osteoarthritis). We also detected LOF burdens for different genes at the same locus (PMVK and THBS3). For the same phenotype, the effect sizes in the LOF burden analysis are consistently larger compared with those identified in the common variant analysis, except for VIT, for which they are the same.
We identify eight interconnected biological pathways that are enriched for effector genes, the majority of which are newly reported here (Table 1, Methods, Supplementary Note, Supplementary Tables 13 and 26–29 and Extended Data Fig. 3; a detailed description of these pathways and the role of the effector genes is provided in the Supplementary Note). We find that the biological processes with the highest number of effector genes, such as ECM and WNT signalling, show higher levels of osteoarthritis heritability explained (Supplementary Fig. 9).
The retinoic acid signalling pathway (Extended Data Fig. 4) is associated with the highest-scoring effector gene, ALDH1A2. ALDH1A2 catalyses the synthesis of all-trans retinoic acid (ATRA), which then interacts with retinoic acid and retinoid acid receptors, regulating the expression of multiple genes with fundamental roles in skeletal patterning and differentiation14,15, as well as organ and limb development16,17. CYP26B1 is involved in the degradation of ATRA, thereby controlling its availability. The balance of synthesis and degradation of ATRA is important for receptor interactions, and depletion or excess of ATRA can result in developmental abnormalities18.
TGFβ signalling (Extended Data Fig. 5) is intricately involved in the pathogenesis of osteoarthritis through its effects on chondrocyte and osteoblast differentiation, skeletal development, cartilage and bone formation, inflammation, ECM remodelling, osteophyte and synovial tissue changes, and interactions with other signalling pathways, such as BMP. The identified effector genes traverse all aspects of TFGβ signalling (Extended Data Fig. 5). We find that TGFB1 and SMAD6 demonstrate allelic imbalance in subchondral bone (Methods, Supplementary Table 27 and Supplementary Fig. 10) and that the osteoarthritis risk allele of rs146652543 is associated with decreased expression of TGFB1. We also identify decreased protein abundance of TGFβ1 in degraded compared with intact osteoarthritis cartilage (Supplementary Table 13). The hip osteoarthritis risk-increasing allele of rs2469081 is associated with decreased expression of SMAD6, a newly identified signal. Furin plasma protein quantitative trait loci (pQTLs) colocalize with osteoarthritis signals on chromosome 15 (rs1894401) (Methods and Supplementary Table 28).
BMP signalling has an important role in many organs and tissues during early embryogenesis (dorsoventral and anteroposterior axis formation), and in postnatal homeostasis. The role of BMP signalling in skeletal development and maintenance is well established, with a lack or excess of BMP signalling giving rise to skeletal abnormalities. Mutations and/or deletion of the effector genes BMP2, BMP6, BMPR1B, GDF6 and GDF5 have been associated with brachydactyly (BMP2, BMPR1B and GDF5)19,20,21, joint deformities and osteoarthritis (GDF5)22, reduction in long bone size (BMP6)23, joint defects (GDF5 and GDF6) and severe chondrodysplasia (BMP2)24. The mechanisms of involvement of BMP signalling with osteoarthritis pathology are complex, ranging from embryonic and developmental changes to those that occur throughout life, such as cartilage homeostasis, osteophyte formation and subchondral bone changes.
WNT signalling has an important function in bone and cartilage metabolism and a well-established role in osteoarthritis25. Two of the effector genes involved in this pathway are WNT family members (WNT3 and WNT5a), both newly reported here, and the remaining genes are involved in modulation of the WNT signalling pathway. WNT signalling has an essential role in embryonic development and homeostasis of bone and cartilage. Dysregulated WNT signalling can contribute to various aspects of osteoarthritis pathology, including cartilage degradation, subchondral bone changes, synovial inflammation and osteophyte formation. We find that the hip osteoarthritis risk allele of rs77601616 is associated with increased expression of SFRP4, located at a locus newly discovered here, in subchondral bone (Methods, Supplementary Table 27 and Supplementary Fig. 10).
Members of the fibroblast growth factor (FGF) pathway have been implicated in the pathogenesis of osteoarthritis through skeletal development, bone and cartilage homeostasis, and also through inflammation and angiogenesis. Five of the effector genes involved in FGF signalling are key FGF pathway members (FGF1, FGF18, FGFR3, FGFR4 and FGFRL1). FGFs have an important role in tissue regeneration and repair and are integral to cell differentiation, proliferation, apoptosis, metabolism, morphogenesis and tissue healing. Two FGF-related pathways involve a further 18 effector genes: FGFR3 signalling in chondrocyte proliferation and terminal differentiation (10 effector genes), and osteoarthritic chondrocyte hypertrophy (16 effector genes) (Supplementary Table 29 and Supplementary Fig. 11). Mutations in FGFR3 are known to give rise to achondroplasias26. Osteoarthritic chondrocyte hypertrophy is associated with dysregulation of FGF, hypoxia and angiogenesis27.
Among the 61 effector genes associated with ECM assembly and organization, 14 are collagens, 3 are proteoglycans, 12 are glycoproteins, 6 are ECM secreted factors, 7 are ECM regulators and 1 is an ECM-affiliated protein. The majority of the ECM in healthy articular cartilage is composed of aggrecan, encoded by ACAN, and collagen type II, encoded by COL2A1, both newly reported effector genes. Mutations in both COL2A1 and ACAN give rise to types of spondyloepiphyseal dysplasia characterized by premature osteoarthritis28. During osteoarthritis progression, the balance between the aggrecan content (which provides the ability to withstand compression and absorb shocks) and collagen content (which provides tensile strength) is critical. Changes in ECM content can give rise to reduced mechanical strength, lack of elasticity and increased susceptibility to damage. We find further support for the involvement of COL2A1 for the association signal at rs11168351, which colocalizes with COL2A1 plasma pQTLs (Supplementary Table 28). The pericellular matrix, which surrounds the chondrocyte and modulates the environment, is enriched for collagen type VI (COL6) and perlecan (HSPG2). COL6 is encoded by six genes, two of which are effector genes (COL6A1 and COL6A2). Mutations in COL6A1/2 are associated with various myopathies29. Mutations in HSPG2, which is also an effector gene, give rise to Schwartz–Jampel syndrome type 1, characterized by myotonia and chondrodysplasia30. Two genes involved in the ECM also harbour LOF burdens (COLGALT2 and COL27A1). The LOF + MIS burdens in COLGALT2 are protective against osteoarthritis (Supplementary Table 23). COLGALT2 encodes an enzyme that is involved in the post-translational glycosylation of collagens and proteins containing collagen domains. Differential allelic expression imbalance between intact and degraded cartilage has shown that lower expression of COLGALT2 is protective for osteoarthritis31. In osteoarthritic cartilage, the risk allele of rs11583641 was associated with increased expression of COLGALT2 mediated through decreased methylation32. Mechanistically, over-glycosylation may result in weakened integrity of collagen fibrils and decreased resilience of the cartilage. The risk of disease was increased for LOF variants in COL27A1, which is a fibril-forming collagen with a role in the transition of cartilage to bone during skeletogenesis. COL27A1 has been shown to be regulated by SOX9 (an effector gene). Mutations in COL27A1 are associated with Steel syndrome, characterized by short stature, hip dislocation and scoliosis33,34.
The circadian rhythm has not been genomically linked with osteoarthritis, although a few studies have established a role for circadian clocks in articular cartilage in regulating pathways related to tissue ageing, degeneration and osteoarthritis. It has also been demonstrated that chronic circadian misalignment may accelerate tissue ageing and ECM degradation. Furthermore, changes in tissue stiffness, for example during ageing, can impair circadian clock function35,36,37. A subpopulation of chondrocytes has also been shown to have increased expression of circadian-related genes (PER1 and SIRT1)38. Disruptions to circadian rhythms may affect the ability of bone and joint tissues to repair and regenerate. Morning joint stiffness can occur due to circadian variations, and age-related changes in sleeping patterns can decrease the amplitude of circadian rhythms. Circadian rhythms can also influence pain perception and sensitivity39, and the absorption, distribution and metabolism of drugs. Circadian-related pain perception has been observed in individuals with osteoarthritis of the knee and hand40,41. Effector genes implicated in this biological process are core circadian clock components (CLOCK, ARNTL and NR1D1), involved in clock entrainment, orchestration, sleeping patterns, transcription of clock genes, circadian oscillations and/or clock-controlled autophagy in bone metabolism (Extended Data Fig. 6 and Supplementary Fig. 12). We find that GFPT1, linked with clock entrainment, demonstrates allelic imbalance in subchondral bone; and that the hip osteoarthritis risk allele of rs6546511 is associated with increased GFPT1 expression (Methods, Supplementary Table 27 and Supplementar