Deciphering osteoarthritis genetics across 826,690 individuals from 9 populations

1Osteoarthritis is one of the leading causes of disability and pain worldwide, with over 300 million people affected. Currently no curative treatments are available, and management strategies focus on symptom alleviation through pain relief and arthroplasty. A detailed understanding of disease etiopathology and new drug targets are therefore urgently needed.

2Osteoarthritis is a complex degenerative disease of the whole joint, characterized by cartilage degeneration, subchondral bone thickening, osteophyte formation, synovial inflammation, and structural alterations of the joint capsule, ligaments, and associated muscles. Recently, advances were made in elucidating the genetic background of osteoarthritis, using genome-wide association studies (GWAS), with 96 statistically independent risk variants reported to date. These variants only explain a small proportion of the phenotypic variance and are mainly associated with osteoarthritis affecting the knee and hip joints.

3Osteoarthritis can affect every synovial joint and an increase in body mass index (BMI) is associated with risk of disease. A better understanding of the genetic differences between weight bearing (knee, hip, and spine) and non-weight bearing joints (hand, finger, and thumb) is needed to help disentangle the metabolic and biomechanical effects contributing to disease development. Here, we conducted a GWAS meta-analysis across knee, hip, finger, thumb, and spine osteoarthritis phenotypes in 826,690 individuals of European and East Asian descent. We integrated functional genomics analyses from disease-relevant tissue, including gene expression, protein abundance and genome-wide methylation, mouse knockout model and monogenic human disease phenotyping data, and complementary computational fine-mapping, colocalization, and causal inference approaches to identify likely effector genes and facilitate much-needed translation into therapies by enhancing our understanding of disease etiopathology.

4We performed GWAS meta-analyses for osteoarthritis across 13 international cohorts stemming from 9 populations (Table S1), in up to 826,690 individuals (177,517 osteoarthritis patients). This is a substantial (2.3-fold) increase of osteoarthritis patient numbers compared to the largest osteoarthritis GWAS to date. Two of the cohorts are of East Asian and 11 of the cohorts are of European descent. We defined 11 phenotypes encompassing all major sites for osteoarthritis (Figure 1; Table S1; STAR Methods). We found 11,897 genome-wide significantly associated single nucleotide variants (SNVs) using a threshold of p < 1.3 × 10−8, to account for the effective number of independent tests. We applied conditional analyses within phenotype and identified 223 independent associations, some of which overlap across phenotypes (Figure 1; Table 1). Eighty-four variants have not been associated with osteoarthritis before. We investigated the previously reported osteoarthritis-loci and found that 87 out of 96 replicated in the same direction at nominal significance (Table S2).

5We used conditional analyses to identify associations that do not overlap across disease phenotype definitions. We identified 100 unique and independent variant associations, 60 of which were associated with more than one osteoarthritis phenotype. Fifty-two variants have not been associated with any osteoarthritis phenotype before (Tables 2 and S3). For each of the 100 association signals, we defined the lead SNV as the risk variant with the strongest statistical evidence for association. Six lead SNVs are coding (all missense), 59 reside within a gene transcript, and 35 are intergenic.

6Here, we report signals for spine (n = 1) and thumb (n = 2) osteoarthritis and increase the number of risk SNVs for hand (5 new, 3 previously reported) and finger (3 new, 2 previously reported) osteoarthritis, phenotypes that had not been studied at scale before (Tables 1, 2, and S3). Of the 100 SNVs, 90 are common (minor allele frequency [MAF] ≥5%) and 4 are low-frequency variants (MAF <5% and ≥0.5%). We detected 6 rare variant associations (MAF 0.03%–0.11%) with large effect sizes (odds ratio [OR] range = 3.03–9.52) (Table 2); 1 variant association was previously reported and 5 variant associations are new findings. All of the new rare variant associations are primarily driven by a large extended family in Iceland.

7Signals from 4 osteoarthritis phenotypes (spine, knee, knee and/or hip, and osteoarthritis at any site) included individuals of non-European ancestry (between 0.9%–2.8% of cases were of East Asian descent). Even though sample sizes in the East Asian cohorts are small, we observed that 62% of the signals have supportive evidence in East Asian ancestry-only analysis, with the same direction of effect, and 20% of these signals are also nominally significant (binomial test p = 2.27 × 10−5, 95% confidence interval [CI] = 7%–100%) (STAR Methods).

8We investigated the predictive power of polygenic risk scores (PRS) and found significantly higher odds of developing disease in individuals at the higher decile of the PRS distribution for several osteoarthritis phenotypes (Table S4; STAR Methods).

9To investigate the presence of osteoarthritis signals specific to males only, females only, or with effects of opposite direction in men and women, we performed a sex-differentiated test of association and a test of heterogeneity in allelic effects. We identified 3 new female-specific independent SNVs, two of which showed significant (Phet-diff <0.016) differences in effect size between sexes (Tables 2 and S5). rs116112221 (Psex-diff = 3.20 × 10−9, Phet-diff = 4.09 × 10−4; female OR = 1.95, 95% CI = 1.58–2.41, P-female = 4.61 × 10−10; male OR = 1.06, 95% CI = 0.82–1.38, P-male = 0.64) is significant in the female-only total hip replacement phenotype and is located in a region containing long intergenic non coding RNAs with the closest protein coding gene being FANCL. FANCL mutations are potentially causative for premature ovary insufficiency in humans, a condition that leads to early menopause, which has been suggested to be linked to increased prevalence of osteoarthritis, although definitive evidence for this hypothesis is still lacking. Preclinical and clinical studies indicate that selective estrogen receptor modulators (SERMs) treatment has consistently positive effects on osteoarthritis, especially for postmenopausal patients with early-stage or osteoporotic osteoarthritis.

10We further identified a signal associated with total hip replacement with opposite direction of effects between men and women, rs10282983 (Psex-diff = 4.93 × 10−16, Phet-diff = 7.66 × 10−14; female OR = 1.15, 95% CI = 1.11–1.19, P-female = 2.21 × 10−14; male OR = 0.92, 95% CI = 0.88–0.96, P-male = 5.16 × 10−4). rs10282983 resides in an intron of C8orf34, which has been associated with waist-to-hip ratio and heel bone mineral density, both risk factors for osteoarthritis. rs10453201 is significantly associated with female osteoarthritis at any site (Psex-diff = 5.67 × 10−9, Phet-diff = 0.049; female OR = 1.05, 95% CI = 1.03–1.06, P-female = 1.05 × 10−8; male OR = 1.02, 95% CI = 1.003–1.04, P-male = 0.02) and is located 5′ of UBAP2, which has been associated with Parkinson’s disease, type 2 diabetes, BMI, and heel bone mineral density in humans.

11Genome-wide meta-analysis identified a new risk variant for early osteoarthritis with large effect size and low allele frequency (rs148693048; effect allele frequency = 0.12%, p = 3.37 × 10−8, OR = 6.26, 95% CI = 3.26–12.00) (Tables 2 and S3). The variant is nominally significantly associated in all contributing studies and with the same direction of effect. rs148693048 has not been associated with osteoarthritis before. Two protein-coding genes in the vicinity show significantly different expression levels in intact compared to degraded cartilage (NEFM and DOCK5). NEFM (neurofilament medium) is relevant to the elongation of neuronal structures, and the expressed protein is commonly used as a biomarker of neuronal damage. The guanine nucleotide exchange activity of DOCK5 (dedicator of cytokinesis 5) has been identified as a regulator of osteoclast function, playing an essential role in bone resorption. Pharmacological inhibition of its activity prevents osteolysis, while preserving bone formation in both humans and mice. Intronic variation in DOCK5 also shows association (p < 5.0 × 10−8) with other bone phenotypes, such as heel bone mineral density and adolescent idiopathic scoliosis.

12We observed that some variants demonstrate a joint-specific effect. We found that the majority of SNVs (60 out of the 100) were genome-wide significantly associated with more than one osteoarthritis phenotype (Figure 2). Forty of the identified SNVs show genome-wide significant associations with weight bearing joints only and 4 SNVs show genome-wide significant associations with non-weight bearing joints only (Figure 2; Table S3). We have over 80% power to detect all 4 non-weight bearing specific variants in the weight bearing joint analyses (at genome-wide significance). Further, we have over 80% power to detect 22 of the 40 weight bearing joint-specific effects in non-weight bearing joint analyses (hand osteoarthritis). Although several core pathways are known to underpin osteoarthritis pathology, regardless of joint site affected, no common genetic osteoarthritis SNVs have been found previously, with the exception of the GDF5 locus. Here, we have identified 42 SNVs with strong association across both weight bearing and non-weight bearing joints. Several of these SNVs, rs3771501 (TGFA), rs3993110 (TEAD1/DKK3), rs72979233 (CHRDL2), and rs7967762 (PFKM/WNT10B) (Figures 2B and 2D), are associated with multiple osteoarthritis joint sites. These signals likely represent a common underlying mechanism in osteoarthritis pathology. They have been shown to play a role in the transforming growth factor β (TGF-β)/bone morphogenetic protein (BMP), Wnt/β-catenin signaling pathways, the functional interaction of which has been implicated in the pathogenesis of osteoarthritis. These signaling pathways could be prime candidates for drug development.

13Additional insights may also be gleaned from the comparison of association signals across osteoarthritis phenotypes. Most of the SNVs associated with knee, hip, and knee and/or hip osteoarthritis have a larger effect size on the respective joint replacement-defined phenotypes, all of which are notably of smaller sample size. This could be driven by homogeneity of phenotype definition (Table S1) or can represent a biological and functional relevance, indicating that these loci might play more important roles in receiving a joint replacement (i.e., pain and inflammation) than in osteoarthritis pathology itself. For example, rs76340814 (PTCH1) and rs28929474 (missense variant in SERPINA1) have stronger associations and larger effect sizes with total hip replacement (THR), total knee replacement (TKR), and total joint replacement (TJR), than with hip or knee osteoarthritis (Figure 2A). Indeed, PTCH1 is thought to function in neurogenic and brain development, and SERPINA1 is thought to function in inflammation. Studies in rat osteoarthritis models have shown that early treatment with alpha-1-antiproteinase, encoded by serpina1, blocked the proteolytic activity of neutrophil elastase and caused lasting improvements in joint inflammation, pain, and saphenous nerve damage.

14We found osteoarthritis subtypes to share substantial genetic components, albeit with a wide range (Figure 2E; Table S6).

15We investigated if osteoarthritis genetic components are shared with other traits and found significant correlation with anthropometric traits (BMI, obesity, weight, and fat mass), type 2 diabetes, education, depressive symptoms, smoking behavior, bone mineral density, reproductive phenotypes and intelligence as previously reported, and several pain phenotypes (Table S6).

16Pain is the most disabling symptom experienced by osteoarthritis patients and is one of the main reasons to proceed to physician consultation and total joint replacement. The etiology of pain in osteoarthritis is multifactorial including significant soft tissue inflammation, the sensitization of pain pathways involving the joint nociceptors, the nociceptive processing in the CNS, and neuropathic pain components in osteoarthritis models. Although a main symptom, no genetic determinants of osteoarthritis pain have been discovered before. We found high correlation between osteoarthritis and sciatica, fibromyalgia, headaches, and other back pain phenotypes, where the highest correlation is with spine osteoarthritis (genetic correlation [rg] = 0.61, 0.87, 0.39, and 0.79, respectively). SOX5, one of the new signals, has been previously reported to be upregulated in human osteoarthritis cartilage and has been associated with back pain and with lumbar intervertebral disc degeneration. These findings are supported by animal model data, in which inactivation of SOX5 leads to defects in skeletogenesis such as in cartilage development, the notochord, and intervertebral discs in mice. We also observed strong correlation between osteoarthritis and pain phenotypes in the LD-Hub database (all derived from the UK Biobank dataset), in particular between spine osteoarthritis and dorsalgia (rg = 0.87), leg pain on walking (rg = 0.82), knee pain (rg = 0.63), hip pain (rg = 0.76), back pain (rg = 0.75), and neck/shoulder pain (rg = 0.67) (Table S6). Thus, our data suggest that a proportion of the identified signals are also associated with osteoarthritis pain.

17We employed complementary computational approaches (STAR Methods) to fine-map the GWAS signals to a small set of likely causal variants, identify relevant tissues based on signal enrichment (Figure S1), and provide mechanistic insights based on expression quantitative trait locus (eQTL) colocalization and causal inference analysis (Table S7). Twelve signals were fine-mapped to variant sets contained entirely within the transcript of a single gene with >95% posterior probability (PP), although we note that this does not provide conclusive evidence for the effector gene. Of note, ALDH1A2, which fine-maps to 6 intronic variants with 99% PP, is currently the target of approved drugs in use for other indications, providing a potential opportunity for drug repositioning (Table S8).

18For 6 SNVs (3 new and 3 known), a single variant could be postulated as causal with >95% PP (Table S8).

19We assessed if any of the genes residing within 1 Mb of the osteoarthritis-associated lead variants showed differential gene expression and protein abundance in primary osteoarthritis-affected tissue in chondrocytes extracted from osteoarthritis patients undergoing joint replacement surgery. Similarly, we compared gene expression of subchondral bone tissue underneath the intact and degraded cartilage tissue (Tables S9 and S10). By combining results from the complementary functional genomics and computational approaches (outlined above), we identified 637 genes with at least one line of evidence pointing to a putative effector gene (Table S10). For these 637 genes, we combined supportive information from the fine-mapping, eQTL colocalization analyses, animal model data, human musculoskeletal and neuronal phenotype data, functional genomics, and causal inference analysis and identified 77 genes that have at least 3 different lines of evidence in support of their role as an effector gene (Tables 3 and S10). Of these 77 genes, 4 are supported by missense lead variants (rs2276749 in VGLL4, rs3740129 in CHST3, rs143083812 in SMO, and rs4252548 in IL11). Forty eight provide strong additional evidence for the likely effector gene at previously reported osteoarthritis-associated SNVs (Table 3) and 30 reside in newly associated signals.

20CHST3, SMAD3, and GDF5 accrued the highest levels of confidence, each with 6 different lines of evidence in support of their involvement in osteoarthritis. CHST3 (carbohydrate sulfotransferase 3) represents a newly identified signal and encodes chondroitin sulfate, the major proteoglycan present in cartilage. Mutations in CHST3 have been previously associated with short stature, congenital joint dislocations, clubfoot, Larsen syndrome, and elbow joint dysplasia. CHST3 has also been shown to be associated with lumbar disc degeneration.

21To glean further insight into the biological role of the high-confidence effector genes in disease processes, we integrated additional information based on the analysis of endophenotypes more closely related to the underlying biology, monogenic and rare human disease data, phenome-wide analyses, and additional functional genomics data (Tables S11 and S12; STAR Methods). By synthesizing all lines of evidence, we found that the assignment of several of the 77 high-confidence effector genes into likely mechanisms through which they exert their effect traverses multiple biological processes (Figure 3A). Here, we primarily focus on the newly associated genes that are reported in this work. These represent high-value candidates for further mechanistic and clinical investigation.

22The majority of high-confidence effector genes are associated with skeletal development (63 in total, 21 genes associated with newly reported signals) and joint degradation (50 in total, 18 genes associated with newly reported signals; 13 genes in common between the skeletal development and joint degradation categories) (Figure 3A). Three effector genes arising from new genetic signals encode structural proteins: CHST3, COL2A1, and FBN2. Collagen type II alpha 1 chain (COL2A1) codes for an essential structural component of cartilage and is important for joint formation and bone growth (Figure 3B). A wide spectrum of diseases is associated with COL2A1, including cartilage and bone abnormalities, such as spondyloepimetaphyseal dysplasia, Kniest dysplasia, and early onset osteoarthritis. Fibrillin 2 (FBN2) encodes a glycoprotein that forms microfibrils in the extracellular matrix and has a major role during early morphogenesis. Fibrillins potently regulate pathways of the immune response, inflammation, and tissue homeostasis, are important in bone remodeling, and regulate local availability of BMP and TGF-β (Figure 3B). Mutations in FBN2 cause contractual arachnodactyly.

23Several genes are connected with signaling pathways. Vestigial like family member 4 (VGLL4) functions via interacting with TEA domain (TEAD) transcription factors. Notably, we identified another new THR and hand osteoarthritis-associated signal located in such a transcription factor, the TEAD1 gene, indicating a common molecular pathway underlying both signals (Figure 3B). TEAD1 functions in the Hippo signaling pathway and is transcriptionally regulated by the YAP1 and TAZ protooncogene proteins, which are involved in mechanosensing and mechanotransduction. Mechanoadaptation of articular cartilage is an important factor in osteoarthritis. Downregulation of VGLL4 is linked to the upregulation of Wnt/β-catenin pathway target genes.

24Wnt family member 1 (WNT1) and wnt family member 10B (WNT10B) are involved in the Wnt signaling pathway, which has an established role in osteoarthritis pathogenesis. Mutations in WNT10B have been linked to limb defects and dental abnormalities, and mutations in WNT1 are associated with osteogenesis imperfecta. Insulin like growth factor 1 receptor (IGF1R) has tyrosine kinase activity, mediates the action of insulin-like growth factor, and regulates cartilage mineralization.

25Nitric oxide synthase 3 (NOS3) encodes the vascular endothelium isoform of nitric oxide synthase (eNOS). NOS3 is associated with sporadic limb defects in mice and has been implicated in bone remodeling in rats. LIM homeobox transcription factor 1 beta (LMX1B) is a transcription factor. Mutations in LMX1B cause a rare autosomal dominant disorder characterized by dystrophic nails, hypoplastic or absent patellae, and dysplasia of the elbows and iliac horn.

26Patched 1 (PTCH1) codes for a receptor for Hh ligands and regulates the activity of smoothened, frizzled class receptor (SMO,another effector gene associated with a known lead SNV). When bound, PTCH1 relinquishes its inhibitory effect on SMO and activates the Hh signaling cascade, which plays an important role in controlling the proliferation of chondrocytes and also in stimulating osteogenesis during endochondral bone formation and longitudinal growth.

27Several further newly identified high-confidence effector genes have a neuronal connection (Figure 3A). Augurin, the protein encoded by C2orf40 (also called ECRG4), is involved in CNS development in animal models and shows association with neuropathologic features of Alzheimer’s disease and related dementias in humans. SNVs in the vicinity of TSEN15 have been robustly associated with anthropometric traits that have epidemiological links to osteoarthritis, such as height, body fat distribution, and waist circumference adjusted for BMI. Cut like homeobox 1 (CUX1) is a transcription factor involved in brain neuronal differentiation and synaptogenesis. Cux1 expression was observed at chondrogenic interzones during limb development, suggesting also a regulatory role in joint formation.

28The TRIO and f-actin binding protein (TRIOBP) gene encodes multiple protein isoforms via 2 promoters. TRIOBP-1 is ubiquitously expressed and interacts with TRIO and f-actin binding protein that together play crucial roles in neuronal morphogenesis and controlling actin cytoskeleton organization, cell motility, and cell growth.

29Myotubularin related protein 2 (MTMR2) has an important role in membrane targeting, vesicular trafficking, and regulation of signal transduction pathways. Mutations in MTMR2 cause Charcot-Marie-Tooth disease type 4B, which features a generalized loss of large myelinated nerve fibers and focally folded myelin sheaths giving rise to inadequate nerve signaling to muscles, resulting in muscle weakness and atrophy. The ubiquitously expressed protein encoded by CREB-binding protein (CREBBP) plays a critical role during development in particular with brain size regulation, correct neural cell differentiation, and neural precursor cell migration, as demonstrated in mouse models.

30Cholinergic receptor muscarinic 2 (CHRM2) is involved in the mediation of cellular responses. Analysis of rat tissues revealed expression in whole brain and in human neuroblastoma cells. Variation in CHRM2 predisposes to various neuropsychiatric diseases, and Alzheimer’s disease. The protein encoded by synaptosome associated protein 47 (SNAP47) is a soluble N-ethylmaleimide-sensitive fusion protein attachment protein receptor (SNARE) protein involved in trafficking and membrane fusion. SNARE-mediated fusion is an essential mechanism that drives the synaptic transmission, neuron development, and growth. SNAP47 plays a role in exocytic mode and neuronal morphogenesis.

31Several of the effector genes have an immune or inflammatory role. For example, the protein encoded by toll like receptor 4 (TLR4) plays a fundamental role in pathogen recognition and activation of the innate immune response. TLR4 is also activated by host-derived molecules generated by damaged tissues related to different musculoskeletal pathologies. This, along with gene expression in chondrocytes, osteoblasts, and synoviocytes, has linked TLR4 to diseases like rheumatoid arthritis, osteoarthritis, and osteoporosis, where modulation or inhibition of TLR4 has been suggested as a treatment. Activation of T cells can lead to osteoclastogenesis and bone resorption by influencing the expression of tumor necrosis factor ligand superfamily member 11 (TNFSF11). TNFSF11 encodes receptor activator of nuclear factor kappa-β ligand (also known as RANKL), a cytokine that has been linked to inflammatory bone remodeling in rheumatoid arthritis, with increased TNFSF11 levels associated with worsening arthritis severity and a well-established role in osteoclastogenesis.

32Nuclear receptor subfamily 3 group C member 1 (NR3C1) encodes the glucocorticoid receptor (GR) which circulates in the cytoplasm and is involved in the inflammatory response. In osteoarthritis, endogenous glucocorticoid signaling in osteoblasts and chondrocytes is detrimental.

33Phosphofructokinase (PFKM) has a role in muscle function. It encodes a muscle isozyme that catalyzes the phosphorylation of fructose-6-phosphate during glycolysis. Mutations in this gene result in Tarui’s disease (glycogen storage disease type 7) that is an autosomal recessive metabolic disorder characterized clinically by exercise intolerance, muscle cramping, exertional myopathy, and compensated hemolysis.

34We examined the druggability status of all 637 genes with at least one piece of supporting evidence from fine-mapping and functional analyses (Table S10; STAR Methods). Of these 637 genes, 205 were present in the druggable genome database, showing a 1.46-fold enrichment of genes with supporting evidence in the database (binomial test p = 2.21 × 10−9) (STAR Methods). From these osteoarthritis druggable target genes, 71 genes reside in tier 1, which incorporates the targets of approved (licensed) drugs and drugs in clinical development (Table S10; STAR Methods). Of the 77 genes with three different lines of evidence supporting causality, 20 are tier 1 candidates (18 of these are present in DrugBank) (Table 4; STAR Methods), of which 7 correspond to new genetic signals discovered in this study (CHST3, VDR, TNFSF11, IGF1R, NR3C1, CHRM2, and NOS3).

35Within tier 1, ten candidates have previously been studied in clinical trials of efficacy or in cohort studies of osteoarthritis (six arising from new signals: PPARD, NR3C1, VDR, MAPK14, IGF1R, and CHST3). The PPARD antagonist sulindac has marketing authorization as a non-steroidal anti-inflammatory drug (NSAID) in osteoarthritis for its prostaglandin synthase activity. The SLC1A1 agonist and neuropathic pain inhibitor pregabalin is commonly prescribed in osteoarthritis. Pregabalin has some supportive clinical trial data for its co-prescription with the NSAID meloxicam in the short-term treatment of pain in knee osteoarthritis. NR3C1 encodes the glucocorticoid receptor, the activation of which has broad anti-inflammatory and immunomodulatory actions with marketing authorization for several agonist molecules. One of these, prednisolone, has long been used as a disease modifying agent in inflammatory arthritis and in the recent Heart Outcomes Prevention Evaluation (HOPE) study was found to be effective in reducing pain and synovitis in hand osteoarthritis. Cathepsin K (encoded by the CTSK gene) is an enzyme that plays a critical role in collagen degradation within osteoclasts, and MIV-711 is a selective cathepsin K inhibitor that has recently been shown in a phase 2 clinical trial to be effective in reducing structural damage in patients with knee osteoarthritis. VDR encodes the vitamin D receptor, the activation of which is a major regulator of calcium metabolism. The results of clinical trials of vitamin D supplementation on symptoms and structural damage in knee osteoarthritis have been mixed but may suggest a small benefit in patients with vitamin D deficiency. EGLN2 encodes Egl nine homolog 2, a prolyl hydroxylase that mediates hydroxylation of proline and thus contributes to collagen and proteoglycan synthesis. Supplementation of its agonist, ascorbic acid (vitamin C), has been associated with joint health in observational cohorts, although with mixed effects. Deficiency of the HCAR2 agonist niacin (vitamin B3) was associated with knee osteoarthritis progression in the Japanese ROAD cohort. The MAPK14 antagonist PH-797804 has been studied in a phase 2 clinical trial to examine the pain relief of PH-797804 alone or with naproxen in subjects with osteoarthritis of the knee (NCT01102660), although we are not aware of any trial results reporting in PubMed or on ClinicalTrials.gov. Finally, the carbohydrate sulfotransferase 3 agonist thalidomide has been shown to attenuate early osteoarthritis development in a mouse medial meniscus destabilization model through a mechanism involving the downregulation of vascular endothelial growth factor (VEGF) expression.

36All of the 45 further tier 1 druggable targets have market authorization or are in clinical development for other indications (Table 4). Ten of these are high-confidence effector genes and 16 arise from new genetic signals. The functional and epidemiological evidence of their roles in clinical osteoarthritis presented here provides support for early repurposing investigation. One antibody small molecule, fostamatinib, appears multiple times in Table 4 as a tyrosine kinase inhibitor that targets AAK1, EPHA5, GAK, GSK3A, MAP2K1, MAP2K6, PAK1, and PRKCD, and has marketing approval as a biologic disease-modifying anti-rheumatic drug (DMARD). The JAK2 antibody baricitinib and the TYK2 antibody tofacitinib are both marketed as biologic DMARDs, and the MAP2K1 antibody binimetinib is currently in phase 3 clinical trials as a biologic DMARD. Each of these drugs therefore present a clinical opportunity and putative mechanism for repurposing studies in osteoarthritis. Of the remaining tier 1 and tier 2 druggable targets, the potentially actionable molecules are at an earlier stage of development and present a more distant repurposing opportunity.

37Our findings have generated further knowledge on the differences between weight bearing and non-weight bearing joints and point to mechanisms that are common to disease development at any joint, and joint-type-specific. Indeed, bone and cartilage development pathways were enriched in signals traversing weight bearing and non-weight bearing joints, identifying joint development as a common mechanism for any form of osteoarthritis (Table S13).

38We have been able to establish molecular links between the disease and its main symptom, pain. We demonstrate genetic correlation between osteoarthritis and pain-related phenotypes and identify signal enrichment in neurological pathways (Table S13). Furthermore, several of the high-confidence effector genes have a role in neuropathology. The majority of osteoarthritis cases in this study were defined as total joint replacement and/or self-reported osteoarthritis, and both of these disease phenotypes are highly driven by pain. Identification of these genes can also have implications for further joint pain-related disorders, for which insights have been limited to date.

39A large number of the high-confidence effector genes converge on the endochondral pathway, playing an essential role in homeostasis of the chondrocyte (Figure 3) and osteophytosis. Several of the identified genes are important in TGF-β signaling and function. The newly identified fibrillin 2 (FNB2) signal, together with LTBP1 and LTBP3, regulate the availability of active TGFB1. TGFB1 is the major form of TGFB in cartilage and can activate a cascade of downstream genes through SMAD3-signaling, including ECM-genes which have been identified in our current study, such as carbohydrate sulfotransferase 3 (CHST3).

40Our data provide evidence for the FGF-signaling cascade (FGFR3, FGF18, and PIK3R1) being causally involved in osteoarthritis (Figure 3). FGF18 is currently being tested in clinical trials for its effectiveness in osteoarthritis. The newly identified molecular player phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1) encodes the p85a, p55a, and p50a regulatory subunits of class IA phosphatidylinositol 3 kinases (PI3Ks), which are known to play a key role in the metabolic actions of insulin and are required for adipogenesis. Mutations in PIK3R1 cause agammaglobulinemia 7, immunodeficiency 36, and SHORT syndrome, which is characterized by short stature, hyperextensibility of joints, ocular depression, Rieger anomaly, and teething delay. The balance between chondrocyte proliferation, differentiation, and hypertrophic conversion is controlled by crosstalk between several signaling pathways, of which we find causal evidence here: PTHLH and IHH-signaling (SMO and PTCH1) antagonize signaling through FGFR3. In addition, we identify two independent genetic variants implicating noggin (NOG) as an osteoarthritis effector gene. Noggin binds to TGFB, BMPs, and GDF5 and thereby prevents binding to the cognate receptor. Mutations in NOG cause a whole range of bone and cartilage phenotype depending on the severity of the mutation.

41Several of the putative osteoarthritis causal genes are involved in developmental pathways (Figure 3). Skeletal development can be linked to osteoarthritis in several ways. First, skeletal developmental genes are involved in joint (tissue) characteristics before onset of disease such as cartilage thickness (TGFA, FGFR3, RUNX2, and PIK3R1) or joint shape (resulting in different loading of the joint). Second, the skeletal development pathway could be involved in the reaction to damage in the joint. Depending on the specific genetic makeup, each individual reaction to a damaging trigger to the joint could be different, thereby determining the risk of developing osteoarthritis upon trauma or mechanical overload. Pathway analysis performed on the current study signals further corroborated this, because it revealed evidence for enrichment of pathways typically involved in reaction to damage.

42Our data also suggest that subtle changes in pivotal osteochondrogenic pathways lead to an adverse response to joint damage and/or overload. This may catalyze a fibrotic response both in cartilage and in the synovium. We identified tenascin C (TNC) as one of the high-confidence effector genes (Figure 3). TNC is a component of the extracellular matrix and is involved with organ fibrosis, inflammation, and cardiovascular disease. The formation of fibrocartilage and fibrosis in the joint is a major contributor to the degenerative changes in osteoarthritis. Further, elevated levels of TGF-β signaling are associated with the pathological and fibrosis changes. TGF-β is also a potent inducer of epithelial-mesenchymal transition (EMT). EMT, a process whereby fully differentiated epithelial cells undergo transition to a mesenchymal phenotype giving rise to fibroblasts, is a driver of early fibrosis, which is a typical response to injury or pathological changes and inflammation, all common endpoint outcomes in osteoarthritis. The severity of fibrosis contributes to the degree of degenerative changes that lead to pain in osteoarthritis. We have identified significant association of variants in many of the genes involved in the induction (e.g., EMT genes, CUX1 and multiple molecular components of the TGF-β pathway) and progression of fibrosis (ECM genes e.g., TNC, TGF-β signaling FNB2, LTBP1, LTBP3, TGFB1, and SMAD3) (Figure 3B). These findings indicate that combined variation in the regulation of these genes may collectively contribute to the susceptibility and severity of degenerative changes in osteoarthritis.

43Seventy-one of the implicated genes code for molecules that are the targets of approved (licensed) drugs and drugs in clinical development. Our findings substantially strengthen the evidence for these potential therapeutics, provide drug repositioning opportunities, and offer a solid basis on which to develop, or repurpose, such interventions for osteoarthritis.

44Our work provides a robust springboard for follow-up functional and clinical studies. We have demonstrated clear differences between distinct osteoarthritis patient populations, for example based on disease severity, joint site affected, and sex. We enhance our understanding of the genetic etiology of disease, shed biological insights, and provide a stepping stone for translating genetic associations into osteoarthritis drug development, ultimately helping to catalyze an improvement in the lives of patients suffering from osteoarthritis.

45Enhancing population diversity in genetic association studies is important for discovering risk variants, pinpointing likely causal alleles, improving risk prediction, and ensuring the transferability of findings across global populations. In this work, less than 3% of contributing subjects were of non-European ancestry. Going forward, the identification and inclusion of diverse populations in osteoarthritis genetic association studies is urgently needed.

46Disentangling mechanisms that are active at the point of disease initiation versus those activated during the natural history of disease warrant animal model studies in which disease dynamics can be studied in depth. Indeed, investment in mechanistic studies of the newly identified high-value targets will be important next steps. Clinical trials of intervention will be needed to take our findings forward into mechanism and clinical outcome, therefore elucidating how to target the implicated genes and proteins, how downstream events will be affected, and, ultimately, how these interventions will affect disease outcome in the patient.

47Further information requests should be directed to the lead contact, Eleftheria Zeggini (eleftheria.zeggini@helmholtz-muenchen.de).

48This study did not generate new unique reagents.

49Detailed information relating to the human subjects for each cohort, is provided in Table S1. This includes the number of individuals in each cohort, sex information and ethnicity.

50Arthritis Research UK Osteoarthritis Genetics (arcOGEN) is a collection of 7410 unrelated, UK-based individuals of European ancestry with knee and/or hip osteoarthritis from the arcOGEN Consortium. Samples were collected in 2 stages from 10 United Kingdom locations (London, Nottingham, Oxford, Sheffield, Southampton, Edinburgh, Newcastle-Upon-Tyne, Sheffield, Wansbeck, and Worcester). The majority of cases had primary OA requiring joint replacement of the hip or knee while a smaller number were ascertained by radiographic evidence of disease (Kellgren-Lawrence (KL) grade ≥ 2). The exclusion criteria included the need for joint replacement due to fracture, secondary OA of any cause, and developmental, vascular, or infective causes of joint disease.

51UKHLS also known as Understanding Society, is a longitudinal panel survey of 40,000 UK households (England, Scotland, Wales and Northern Ireland) representative of the UK population. Participants are surveyed annually since 2009 and contribute information relating to their socioeconomic circumstances, attitudes, and behaviors via a computer assisted interview. The study includes phenotypical data for a representative sample of participants for a wide range of social and economic indicators as well as a biological sample collection encompassing biometric, physiological, biochemical, and hematological measurements and self-reported medical history and medication use. https://www.understandingsociety.ac.uk/.

52The ARGO-Larissa study was set up to investigate the genetic architecture of knee OA in a Greek population. It included individuals with primary knee OA undergoing total knee arthroplasty. The osteoarthritis participants’ recruitment was conducted in the city of Larissa, central Greece.

53The ARGO study, was set up to investigate the genetic architecture of hip and knee OA in a Greek population. More than 1,500 patients with severe OA, undergoing hip and/or knee total joint replacement were recruited from the cities of Athens and Larissa. The ARGO collection was conducted in three public hospitals (Attikon University General Hospital of Athens, Nea Ionia General Hospital Konstantopouleio, and KAT Hospital) and one private hospital (Lefkos Stavros General Hospital) in the city of Athens, Greece between February of 2015 and March 2017.

54The UK Biobank study is a large population-based prospective study of > 500,000 participants with ages ranging 40–69 years. In total, 503,325 participants who registered in the National Health Service were recruited out of 9.2 million mailed invitations between 2006 and 2010 in 22 assessment centers throughout the UK. Most participants visited the center once, but some individuals visited the center at up to three times. Baseline data were collected using electronic signed consent, a self-completed touch-screen questionnaire, a brief computer-assisted interview, physical and functional measures, and collection of biological samples and genetic data. The UK Biobank genetic data contains genotypes for 488,377 participants. All detailed genotyping, quality control, and imputation procedures are described at the UK Biobank website (https://biobank.ctsu.ox.ac.uk). Briefly, 50,000 samples were genotyped using the UKBiLEVE array and the remaining samples were genotyped using the UK Biobank Axiom array (Affymetrix) for ∼800,000 SNPs. Population structure was captured by principal component analysis on ∼500,000 UK Biobank samples using ∼100,000 SNPs. After sample and SNP quality control (QC) of the directly-typed genotypes, resulting in 670,739 autosomal markers in 487,442 individuals, data were prephased using SHAPEIT3 and imputed using the IMPUTE4 program (https://jmarchini.org/software/). Both analyses were carried out centrally and the full dataset consisted of approximately 96 million variants in 487,411 individuals. https://www.ukbiobank.ac.uk/

55This work was based on the third UK Biobank release, which includes the full set of genotypes imputed on the Haplotype Reference Consortium and the merged UK10K and 1000 Genomes phase 3 reference panels1000 Genomes Consortium. Access to UK Biobank genetic and phenotypic data was given through the UK Biobank Resource under application request.

56Hong Kong Degenerative Disc Disease Population Cohort (HKDDDPC) is a a population-based cohort with subjects openly recruited via newspapers advertisement, posters and e-mails, regardless of their social and economic status. The study call was for any participant who agreed to a study on the lumbosacral spine with MRI, clinical questionnaires and follow-up assessments. Participants with prior surgical treatment of the spine, spinal tumors and fractures, and marked spinal deformities were excluded from the study. Subjects selected were not based on the presence or absence of clinical symptoms. All qualified subjects underwent T1-weighted axial MRI and T2-weighted sagittal MRI of the lumbosacral spine (L1-S1) after informed consent was obtained from participants and ethics was approved by a local institutional board. MRI Protocol: 1.5T or 3T MRI machines were used for axial and sagittal imaging at L1-S1. Subjects were oriented in the supine position. For T1-weighted axial scans, the field of view was 21cmx21cm, slice thickness was 4mm, slice spacing was 0.4mm, and imaging matrix was 218x256. For T2-weighted sagittal scans, the field of view was 28cmx28cm, slice thickness was 5mm, slice spacing was 1mm, and imaging matrix was 448x336. The repetition time for T1- and T2-weighted MRI were 500ms-800ms and 3320ms respectively, and their echo time was 9.5ms and 85ms. According to the pedicle and disc levels, 11 parallel slices were made at each spinal level with reference to the adjacent endplates. Definition of OA on X-rays: For the lumbar lateral radiographs, each disc level (L1-2, L2-3, L3-4, L4-5, L5-S1) was measured for osteophytes and vertebral narrowing. Grade 0 was considered none; grade 1 was mild; grade 2 was moderate; and grade 3 was severe. The L5-S1 disc was narrowed when its height was less than the disc space of L3-4. Diagnosis of OA was based on the criteria as described by. Disc space narrowing was present with grade ≥ 1 and osteophytes with grade ≥ 2. OA had either “narrowing,” “osteophytes” or “both.” Narrowing was considered with grade ≥ 1 narrowing at 2 or more vertebral levels and osteophytes were considered with grade ≥ 2 at 2 or more vertebral levels. With both narrowing and osteophytes, then “both” was considered. Conversion of MRI for OA diagnosis: Conversion of MRI for OA diagnosis was dependent on two MRI ratings: disc bulging and Schneiderman score. Disc bulging was divided into 4 categories: 0 = no disc bulging; 1 = posterior disc bulging (disc displaced beyond a virtual line connecting the posterior edges of two adjacent vertebrae); 2 = disc extrusion (distance between the edge of the protruded disc into the spinal canal was greater than the distance between edges of the base of the disc); 3 = disc sequestration (disc material detached and migrated away from the level of the intervertebral disc). The Schneiderman score was used to describe the disc signal intensity and was evaluated by a 4-point scale: 0 = normal disc height and signal intensity; 1 = speckled pattern or heterogeneous decreased disc signal intensity; 2 = diffuse loss of signal; 3 = signal void.

57Each lumbar intervertebral disc was rated for disc bulging and Schneidermann score. Lumbar spine OA was defined to be present if (1) At least 1 disc with Schneidermann score 3, OR (2) At least 1 disc with Schneidermann score 2 and Bulging score 2, OR (3) At least 2 discs with Schneidermann score 2 AND disc bulging score 1,

58There were 587 of subjects with X-rays as well. We validated our definition based on MRI by testing and confirming its equivalence with X-ray diagnoses using these subjects. DNA samples were genotyped using the Illumina humanOmniZhongHua-8 v1.2 BeadChip. Quality control (QC) of the genotyped data were conducted based on pipeline provided by. Imputation of single nucleotide polymorphisms (SNPs) was performed using reference panels from the Haplotype Reference Consortium (HRC).

59The Nord-Trøndelag Health Study (HUNT) is a large population-based cohort from the county Nord-Trøndelag in Norway. All residents in the county, aged 20 years and older, have been invited to participate. Data was collected through three cross-sectional surveys, HUNT1 (1984-1986), HUNT2 (1995-1997) and HUNT3 (2006-2008), and has been described in detail previously, with the fourth survey recently completed (HUNT4, 2017-2019). DNA from whole blood was collected from HUNT2 and HUNT3, with genotypes available from 71,860 participants. All genotyped participants have signed a written informed consent regarding the use of data from questionnaires, biological samples and linkage to other registries for research purposes. In total, DNA from 71,860 HUNT samples was genotyped using one of three different Illumina HumanCoreExome arrays (HumanCoreExome12 v1.0, HumanCoreExome12 v1.1 and UM HUNT Biobank v1.0). Samples which failed to reach a 99% call rate, had contamination > 2.5% as estimated with BAF Regress, large chromosomal copy number variants, lower call rate of a technical duplicate pair and twins, gonosomal constellations other than XX and XY, or whose inferred sex contradicted the reported gender, were excluded. Samples that passed quality control were analyzed in a second round of genotype calling following the Genome Studio quality control protocol described elsewhere. Genomic position, strand orientation and the reference allele of genotyped variants were determined by aligning their probe sequences against the human genome (Genome Reference Consortium Human genome build 37 and revised Cambridge Reference Sequence of the human mitochondrial DNA; http://genome.ucsc.edu) using BLAT. Variants were excluded if their probe sequences could not be perfectly mapped to the reference genome, cluster separation was < 0.3, Gentrain score was < 0.15, showed deviations from Hardy Weinberg equilibrium in unrelated samples of European ancestry with p value < 0.0001), their call rate was < 99%, or another assay with higher call rate genotyped the same variant. Imputation was performed on the 69,716 samples of recent European ancestry using Minimac3 (v2.0.1, https://genome.sph.umich.edu/wiki/Minimac3) with default settings (2.5 Mb reference based chunking with 500kb windows) and a customized Haplotype Reference consortium release 1.1 (HRC v1.1) for autosomal variants and HRC v1.1 for chromosome X variants. The customized reference panel represented the merged panel of two reciprocally imputed reference panels: (1) 2,201 low-coverage whole-genome sequences samples from the HUNT study and (2) HRC v1.1 with 1,023 HUNT WGS samples removed before merging. We excluded imputed variants with Rsq < 0.3 resulting in over 24.9 million well-imputed variants.

60Geisinger is an integrated health care provider located in central and northeastern Pennsylvania and New Jersey. Geisinger’s electronic health record (EHR) consists of comprehensive longitudinal clinical information including patients’ demographic data, diagnoses (including co-morbidities), lab measurements, prescriptions, procedures, vital signs, and, of relevance for this study, surgical procedure logs. The EHR captures a median of 14 years of health data for patients within the MyCode® Community Health Initiative biorepository. Through the Geisinger-Regeneron DiscovEHR collaboration, whole exome sequence and genome wide genotype data are available from more than 92,000 MyCode® participants to date. These high dimensional clinical data linked to genetic data provide opportunities to conduct precision health research at an unprecedented scale that can lead to significant clinical insights. https://www.geisinger.org/precision-health/mycode

61The details of MyCode Community Health Initiative have been described previously.

62Genotyping was performed in two batches on the Illumina Infinium OmniExpress Exome array and GSA-24v1-0 array for Geisinger 60k and Geisinger 30k cohorts, respectively. The Michigan Imputation Server was used to impute genotypes for both cohorts to HRC.r1-1 EUR reference genome (GRCh37 build) separately. Pre-imputation QC included sample call rate and marker call rate > 90%, HWE p value > 1e-15, MAF > 1%. A/T & G/C SNPs were removed if MAF > 0.4. SNPs with differing alleles, SNPs with > 0.2 allele frequency difference, SNPs not in HRC reference panel were also removed. Variants with imputation info score > 0.3 and MAC > 5 were included in the analyses. We used ICD-code based method to define OA cases and controls. We adopted a linear mixed model built in BOLT_LMM (N > 5000) or GEMMA (N < 5000) for the association tests while accounting for the relatedness and population structure (first 20PCs). PLINK1.9 was used for genetic data quality control and PC calculation.

63The Japanese cohort of knee OA GWAS (disease cohort) consists of 900 cases and 3,400 controls. The cases were all symptomatic OA. They were diagnosed and recruited by expert orthopedic surgeons based on clinical and radiographic examination. All had clinical records for OA and radiographs (standing knee A-P). The controls were obtained from Biobank Japan. The genotyping was done by using Illumina HumanHap550v3 Genotyping BeadChip. After excluding cases with call rate of < 0.98, we applied SNP QC (call rate of ≥ 0.99 in both cases and controls and P value of Hardy-Weinberg equilibrium test of ≥ 1.0 × 10−6 in controls). Finally, 459,393 SNPs on autosomal chromosomes passed the QC filters.

64The deCODE genetics osteoarthritis study is an ongoing population based study in Iceland that was initiated in 1997. The study includes information on all subjects who have undergone total joint replacement in Iceland, and on osteoarthritis status from the Landspitali University Hospital electronic health records. Information on hand osteoarthritis patients is derived from a database of hand osteoarthritis patients that was initiated in 1972. Subjects have given blood or buccal samples to deCODE genetics biobank, which has gathered genotypic and medical data from more than 160,000 volunteer participants. https://www.decode.com/.

65The details of OA definition and analyses have been described previously.

66The Rotterdam study is a large longitudinal population-based cohort study designed to study the risk factors for all major diseases of the elderly. The study started in 1991, has grown to up to 15000 individuals and has detailed phenotyping for cardiovascular, neurodegenerative, endocrine and locomoter diseases and more. For osteoarthritis, longitudinal X-rays on multiple joints and knee MRI’s are available (and scored), as well as information on joint pain. GWAS data are available for all individuals of the Rotterdam Study, as well as additional molecular layers (such as RNA, methylation, microbiome). Genomic studies in the Rotterdam Study are led by the Genetic Laboratory, Department of Internal Medicine of the ErasmusMC in Rotterdam. http://www.epib.nl/research/ergo.htm, http://www.glimdna.org

67Nurses’ Health Study and Nurses’ Health Study II: The Nurses’ Health Studies are among the largest prospective investigations into the risk factors for major chronic diseases in women (https://nurseshealthstudy.org/). The NHS is a prospective cohort study established in 1976. Blood samples were collected from a subset of participants in 1989-90. DNA was collected from cheek cells from another subset of participants in 2001-2004. The NHS II was established in 1989 to study a population younger than the original NHS cohort. Blood samples were collected on a subset of participants in 1996-1999. DNA from cheek cells was collected in 2006 from another subset of participants. Self-reported cases of total hip replacement from the NHS and NHS2 were analyzed for the GO meta-analysis.

68TwinsUK is the UK’s largest adult twin registry and the most clinically detailed in the world. Professor Tim Spector from King’s College London set up the cohort in 1992 to investigate the incidence of osteoporosis and other rheumatologic diseases in several hundred monozygotic (identical) twins. We now have almost 14000 identical and non-identical twins from across the UK, with ages between sixteen and one hundred and our research has expanded to include multiple diseases and conditions. TwinsUK aims to investigate the genetic and environmental basis of a range of complex diseases and conditions. Current research includes the genetics of metabolic syndrome, cardiovascular disease, the musculoskeletal system, aging, sight as well as how the microbiome affects human health. The TwinsUK cohort is now probably the most genotyped and phenotyped in the world. TwinsUK data have enabled multiple collaborations with research groups worldwide and the publication of research papers. https://twinsuk.ac.uk/

69Details on OA definition as in and GWAS, QC imputation as described in and.

70The GARP study is a prospective observational study in patients with familial generalized osteoarthritis, hand osteoarthritis and other osteoarthritis phenotypes. All patients (N = 380) have symptoms and definite radiological signs of osteoarthritis and represent an advanced disease state. https://www.lumc.nl/org/reumatologie/research/artrose/9031609354853/ The Leiden Longevity Study (LLS) consists of 420 Caucasian families with at least two long-lived siblings (men aged 89 years or above; women aged 91 years or above), the middle aged offspring and the partners of this offspring. https://www.leidenlangleven.nl/en/home

71EGCUT has 52000 gene donors, who are all genotyped. In 2019, 100000 new donors will be collected and genotyped. All osteoarthritis cases were selected from Estonian Biobank which is a population-based biobank of the EGCUT. https://genomics.ut.ee/en

72Osteoarthritis cases were chosen from ca 50000 participants of Estonian Biobank by using ICD 10 codes. To specify THR, TKR and TJR cases, the codes from NOMESCO Classification of Surgical Procedures were used in addition. The numbers for cases and controls used are provided in the attached table.

73For genotyping, Illumina Human CoreExome, OmniExpress, 370CNV BeadChip and GSA arrays were used. Quality control included filtering on the basis of sample call rate (< 98%), heterozygosity (> mean ± 3SD), genotype and phenotype sex discordance, cryptic relatedness (IBD > 20%) and outliers from the European descent based on the MDS plot in comparison with HapMap reference samples. SNP quality filtering included call rate (< 99%), MAF (< 1%) and extreme deviation from Hardy–Weinberg equilibrium (p < 1 × 10−4). Imputation was performed using SHAPEIT2 for prephasing, the Estonian-specific reference panel and IMPUTE2 with default parameters. Association testing was carried out with snptest-2.5.2, adjusting for 4 PCs, arrays, current age and sex (when relevant).

74Avon Longitudinal Study of Parents and Children (ALSPAC). Pregnant women resident in Avon, UK with expected dates of delivery 1st April 1991 to 31st December 1992 were invited to take part in the study. The initial number of pregnancies enrolled is 14,541 (for these at least one questionnaire has been returned or a “Children in Focus” clinic had been attended by 19/07/99). Of these initial pregnancies, there was a total of 14,676 fetuses, resulting in 14,062 live births and 13,988 children who were alive at 1 year of age. When the oldest children were approximately 7 years of age, an attempt was made to bolster the initial sample with eligible cases who had failed to join the study originally. As a result, when considering variables collected from the age of seven onward (and potentially abstracted from obstetric notes) there are data available for more than the 14,541 pregnancies mentioned above. The number of new pregnancies not in the initial sample (known as Phase I enrolment) that are currently represented on the built files and reflecting enrolment status at the age of 24 is 913 (456, 262 and 195 recruited during Phases II, III and IV respectively), resulting in an additional 913 children being enrolled. The phases of enrolment are described in more detail in the cohort profile paper and its update (see footnote 4 below). The total sample size for analyses using any data collected after the age of seven is therefore 15,454 pregnancies, resulting in 15,589 fetuses. Of these 14,901 were alive at 1 year of age. A 10% sample of the ALSPAC cohort, known as the Children in Focus (CiF) group, attended clinics at the University of Bristol at various time intervals between 4 to 61 months of age. The CiF group were chosen at random from the last 6 months of ALSPAC births (1432 families attended at least one clinic). Excluded were those mothers who had moved out of the area or were lost to follow-up, and those partaking in another study of infant development in Avon.

75Ethical approval for the study was obtained from the ALSPAC Ethics and Law Committee and the Local Research Ethics Committees. Informed consent for the use of data collected via questionnaires and clinics was obtained from participants following the recommendations of the ALSPAC Ethics and Law Committee at the time.

76The arcOGEN study was ethically approved by appropriate review committees, and the prospective collections were approved by the National Research Ethics Service in the United Kingdom. All subjects in this study provided written, informed consent.

77The UKHLS has been approved by the University of Essex Ethics Committee, and informed consent was obtained from every participant.

78Verbal informed consent was given by all research participants prior to the collection of blood samples for the research. The research participant recruitment, consent process, and study protocol were approved by the Institutional Review Board of the University Hospital of Larissa and conform to the ethical principles set out in the Declaration of Helsinki (1975).

79The ARGO collection was conducted in three public hospitals (Attikon University General Hospital of Athens, Nea Ionia General Hospital Konstantopouleio, and KAT Hospital) and one private hospital (Lefkos Stavros General Hospital) in the city of Athens, Greece between February of 2015 and March 2017. All studies were approved by the relevant hospital Institutional Review Board and conducted in accordance with the principles set out in the Declaration of Helsinki. All patients provided written informed consent prior to participation.

80All participants signed consent to participate in UK Biobank and UK Biobank’s scientific protocol and operational procedures were reviewed and approved by the North West Research Ethics Committee (REC reference number 06/MRE08/65), North West Multicenter Research Ethics Committee (REC reference 11/NW/0382), the National Information Governance Board for Health and Social Care and the Community Health Index Advisory Group.

81Informed consent was obtained from participants and ethics was approved by a local institutional board.

82All genotyped participants have signed a written informed consent regarding the use of data from questionnaires, biological samples and linkage to other registries for research purposes. The current study was approved by the Regional Committee for Medical and Health Research Ethics (REK) 2015/573.

83MyCode Governing Board and an external Ethics Advisory Council approved the study and informed consent was obtained from all subjects as detailed in this reference:.

84The Ethical committee of RIKEN Yokohama Institute approved the study. Informed consent was obtained from all subjects.

85All participants who donated samples gave informed consent and the National Bioethics Committee of Iceland approved the study (VSN 14-148) which was conducted in agreement with conditions issued by the Data Protection Authority of Iceland. Personal identities of the participant’s data and biological samples were encrypted by a third-party system (Identity Protection System), approved and monitored by the Data Protection Authority.

86The Rotterdam Study has been approved by the Medical Ethics Committee of the Erasmus MC (registration number MEC 02.1015) and by the Dutch Ministry of Health, Welfare and Sport (Population Screening Act WBO, license number 1071272-159521-PG). The Rotterdam Study has been entered into the Netherlands National Trial Register (NTR; https://www.trialregister.nl) and into the WHO International Clinical Trials Registry Platform (ICTRP; https://www.who.int/ictrp/network/primary/en/) under shared catalog number NTR6831. All participants provided written informed consent to participate in the study and to have their information obtained from treating physicians. Nurses’ Health Study and Nurses’ Health Study II: This study was approved by the Institutional Review Boards of the Harvard T. H. Chan School of Public Health and Brigham and Women’s Hospital. Informed consent was obtained from all subjects for the collection of biospecimens for genotyping and the use of their genotype and de-identified data for research.

87Ethics approval was obtained from the Guy’s and St. Thomas’ Hospital Ethics Committee. Written informed consent was obtained from every participant. RAAK, GARP, LLS: Ethical approval for all studies was obtained from the medical ethics committee of the Leiden University Medical Center (RAAK: P08.239 and P19.013; GARP: P76.98; LLS P01.113) and informed consent was obtained from all participants.

88Ethical approval for the study was obtained from the ALSPAC Ethics and Law Committee and the Local Research Ethics Committees. Specific approval references for each clinic can be found at http://www.bristol.ac.uk/alspac/researchers/research-ethics/. Informed consent for the use of data collected at clinics was obtained from participants following the recommendations of the ALSPAC Ethics and Law Committee at the time. Consent for biological samples has been collected in accordance with the Human Tissue Act (2004).

89Our study has been reviewed and approved by the local ethics committee on Estonian Bioethics and Human Research, resolution nr 1.1-12/624. Informed consent was obtained from all subjects.

90Genome-wide association analysis for osteoarthritis was performed across 21 cohorts (Table S1), for a total of 826,690 individuals (177,517 osteoarthritis patients). We defined 11 stratified osteoarthritis phenotypes: osteoarthritis at any site, osteoarthritis of the hip and/or knee, knee osteoarthritis, hip osteoarthritis, total joint replacement, total knee replacement, total hip replacement, hand osteoarthritis, finger osteoarthritis, thumb osteoarthritis and spine osteoarthritis (Figure 1; Table S1). Osteoarthritis was defined by either a) self-reported osteoarthritis, b) clinical diagnosed, c) ICD10 codes (Table S1) or d) radiographic as defined by the TREAT-OA consortium, depending on the data available in the cohort (Table S1). Controls were OA-free or population-based with or without ICD code exclusions. An age exclusion if appropriate (preferable of 45 years and older) was applied at the discretion of each cohort. GWAS analysis were performed by each cohort, and adjusted for cohort specific covariates (Table S1).

91For coding SNVs we considered only the following moderate to high impact annotations when weighting genes for prioritisation: transcript_ablation, splice_acceptor_variant, splice_donor_variant, stop_gained, frameshift_variant, stop_lost, start_lost, transcript_amplification, inframe_insertion, inframe_deletion, missense_variant, protein_altering_variant.

92We investigated if any of the genes within 1Mb (upstream and downstream) of the 100 SNVs had a musculoskeletal or neuronal/pain phenotype in mouse knockouts using information from the Mouse Gene expression database (GDX) of the Mouse Genome Informatics (MGI) database. Mouse orthologs of the genes were extracted from Ensembl, using biomart (GRCh37, Version 69). The MGI Batch Query was used to extract all mouse knockout phenotypes from the GDX for all of the investigated genes Using the MGI mouse phenotype ontology the following mouse knockout phenotypes were included for musculoskeletal phenotypes: skeleton phenotype, muscle phenotype and immune system phenotype. For mouse neuronal/pain phenotype the following MGI mouse phenotype ontology was included: nervous system phenotype. We also investigated if any genes had a musculoskeletal phenotype in mouse knockouts using information from https://www.hugedomains.com/domain_profile.cfm?d=boneandcartilage&e=com. Genes which had a mouse musculoskeletal or neuronal phenotype were reported (Table S10). We also investigated human skeletal genetic disorders. We used the Nosology and classification of genetic skeletal disorders to identify genes within 1Mb (upstream or downstream) of the 100 SNVs that had links to human musculoskeletal phenotypes. For human genes linked to monogenetic pain disorders, we downloaded the OMIM database (https://www.omim.org/) and extracted any genes containing phenotypes with the following key words: pain, pain and neuropathy, neuropathy. In addition, we also included genes linked to pain from the curated Human Pain Genetics Database (https://humanpaingenetics.org/hpgdb/) We then examined if any genes within 1Mb (upstream or downstream) of the 100 SNVs were included in those human pain gene lists. Human genes with a pain phenotype or link were reported (Table S10).

93For additional information on the 100 identified osteoarthritis associated SNVs, we examined their association in several osteoarthritis endophenotype and structural phenotype GWAS studies. The osteoarthritis definition and related structural phenotypes were defined based on radiographs in the Rotterdam Study cohorts of the hip, hand, finger, thumb and knee joints (n = 5,634 to 9,276). We have used the following radiographic measurements to create (semi)-quantitative endophenotypes for the hip, knee, hand, finger and thumb joints: Joint Space Narrowing (JSN) (0–3 scoring), Joint Space Width (JSW) (mm), Osteophytes (0–3 scoring), and Kellgren-Lawrence (KL)-score (0–5). Using these measurements we have defined the following structural OA phenotypes: Finger/Hand/Thumb/Knee/Hip JSN sum score, osteophyte sum score, KL sum score and Hip JSW.JSW was assessed at pelvic radiographs in anterior-posterior position and measured in mm, along a radius from the center of the femoral head. Within the Rotterdam Study, a 0.5 mm graduated magnifying glass laid directly over the radiograph was used to measure the joint space width of the hip joints. Acetabular dysplasia was measured using the Center-Edge angle or also known as the Wiberg (CE-angle). The angle was measured using statistical shape model (SSM) software. A continuous phenotype was used for the CE-angle, because of the normal distribution of the measured angles. Since the CE-angle of the right hip and the left hip has a high correlation (Pearson correlation coefficient 0.68), only the CE-angle of the right hip was used in our GWAS. Minimal Joint Space Width (mJSW) GWAS data was taken from. Summary statistics for Bone Size as measured by DXA were taken from.

94To investigate if any of the high-confidence effector genes show a different expression in osteophytes, indicating a potential role in repair mechanisms in response to joint cartilage degeneration, we investigated if they showed significant (0.1% FDR) differential gene expression, methylation or differential protein abundance in osteophytic cartilage compared to low-grade (intact) cartilage in a within-individual matched analysis from 9 individuals who had undergone THR for primary osteoarthritis (Table S12).

95To identify if any of the high confidence effector genes code for proteins that are implicated in disc degeneration, we investigated if the proteins were differentially abundant in recently published spatiotemporal proteomics atlas of human intervertebral discs. Eight effector genes demonstrate differential protein abundance in a comparison between intervertebral discs from a younger (16 year old male) and an older (59 year old male) individual with no reported scoliosis or degeneration (Table S12).

96We scanned the Human Phenotype Ontology (HPO) database, which is currently being developed using the medical literature, DECIPHER (https://www.deciphergenomics.org/about/overview), OMIM (https://omim.org/) and Orphanet (https://www.orpha.net/consor/cgi-bin/index.php?lng=EN), to examine if any of the high confidence effector genes are implicated in phenotypic abnormalities of monogenic and rare human diseases. Fifty-one genes are involved in diseases related to skeletal development, joint degeneration, adipogenesis, muscle function, neuronal function, immune response and inflammation (Figure 3; Table S12).

97GWAS summary statistics from all cohorts were collected and checked to contain all the data needed for the meta-analysis. The quality control (QC) was performed centrally by using EasyQC. Briefly, missing data, mono-allelic SNVs, nonsensical values (p > 1, infinite beta’s etc.) and duplicates were removed from the data. We excluded variants with poor imputation quality (R2 < 0.3), if the effective sample size was < 20 and if the minor allele count was < 6. Allele coding was harmonized across cohorts (A/T/C/G or I/D). Allele frequency was checked against the imputation reference (HRC http://www.haplotype-reference-consortium.org/) or 1000G http://www.internationalgenome.org/) to identify possible allele coding errors. P values were checked to match the corresponding beta values. Cleaned data was used as input for the meta-analysis. Meta-analysis was performed using inverse variance weighting in METAL. Genomic control was performed on all datasets, except those which had already carried out genomic-control adjustments prior to centralized QC and meta-analysis. Genome-wide significance threshold was set at p < 1.3x10−8, corrected for multiple testing (see significance threshold section below). For each phenotype we only considered variants reported in at least 2 cohorts with the same direction of effect with a minimum MAF > = 0.0001 in any contributing cohort. We repeated the same procedure to perform two sensitivity analyses: a) we excluded from the meta-analyses the largest contributing dataset, UK Biobank (https://www.ukbiobank.ac.uk/) and b) East Asian ancestry-only meta-analyses for the 4 osteoarthritis phenotypes (spine, knee, knee and/or hip, and osteoarthritis at any site) that included East Asian cohorts. To summarize the significance of the signals that have supportive evidence in East Asian ancestry-only meta-analysis, we conducted a binomial test (N = 10 SNVs with concordant direction and p ≤ 0.05, N = 77 SNVs tested, and 0.025 is the expected proportion of SNVs at p ≤ 0.05 and with the same direction of effect). As a sensitivity analysis, we excluded the largest contributing dataset in which the majority of previously-reported loci originate (UK Biobank; up to 68,621 osteoarthritis cases and 247,846 controls) from the meta-analyses.

98The meta-analyses and QC steps described above were repeated for males and females separately in a subset of cohorts. We then combined the resulting association summary statistics (male-specific meta-analysis, consisting of up to 56,462 cases and 153,808 controls, and female-specific meta-analysis, consisting of up to 90,838 cases and 192,697 controls) to conduct a sex-differentiated test of association and a test of heterogeneity in allelic effects, as implemented in GWAMA. This method allows for heterogeneity of allelic effects in magnitude and/or direction between males and females and offers substantial gains in power to detect SNV associations. The genome-wide significance threshold was set at p < 1.3x10−8, corrected for multiple testing. Heterogeneity in allelic effect sizes was assessed with Cochran’s Q statistic and the significance threshold was set at p < 0.016, corrected for the 3 independent new signals identified across the 11 osteoarthritis phenotypes (Table S5).

99We carried out a meta-analysis of early osteoarthritis, defined as age at onset younger than 45 years of age, across 3 cohorts (Estonian Biobank, HUNT & UK-Biobank) with age at onset information available. The analysis was conducted on 6,838 early-onset osteoarthritis patients and 41,449 controls in only one of the 11 phenotype definitions used in the main meta-analysis (osteoarthritis at any site, abbreviated as EarlyAllOA). The QC and meta-analysis steps of the main meta-analysis were repeated for early-onset meta-analysis and the genome-wide significance threshold was set at p < 5x10−8.

100The testing of M = 11 osteoarthritis phenotypes in this study needed to be taken into account in the interpretation of genome-wide statistical significance. Applying a Bonferroni correction would be inherently conservative as this method assumes independence among the tests considered. Therefore, we first used LD Score regression method (https://github.com/bulik/ldsc/) with genome-wide meta-analysis summary statistics to estimate the genetic correlation matrix between the 11 osteoarthritis traits (Table S6) and then calculated the effective number of independent traits (M_eff) from the eigenvalues λ_i of the correlation matrix:For the M = 11 osteoarthritis phenotypes in this study, M_eff = 4.6565. The threshold corrected for the effective number of traits to report genome-wide significance is p < 1.3x10−8.

101To define independent signals for each osteoarthritis phenotype (Table S3), we used the clumping function in PLINK 1.9 with the following parameters: (a) significance threshold for index variants: p < 1.3x10−08, (b) LD threshold for clumping: 0.10, and (c) physical distance threshold for clumping: 1Mb (2Mb window around the index variant). LD calculations were based on the full UK Biobank imputed set. To test that the index variants defined by clumping were statistically independent, we performed an approximate stepwise model-selection procedure, as implemented by COJO in GCTA (https://cnsgenomics.com/software/gcta/). A signal in a region was defined as independent if its P value of association in the stepwise regression was less than the adjusted genome-wide significant threshold (p < 1.3x10−8).

102To define independent signals across the 11 osteoarthritis phenotypes (Tables 2 and S3), we performed reciprocal approximate conditional analyses, as implemented by COJO in GCTA, of each independent variant of one osteoarthritis phenotype conditioned on each independent variant of the other osteoarthritis phenotypes within 1-Mb region. Two signals were considered dependent if the P value for either signal conditioned on the other was either ≥ 1x10−7, or attenuated by at least 2 orders of magnitude. Among dependent variants, the one with the lowest P value was classified as independent. Using an approximate conditional and joint multiple-SNP analysis, as implemented by COJO in GCTA, we investigated the statistical independence between index signals per osteoarthritis phenotype and previously reported osteoarthritis variants within a 1-Mb region. In Tables 1, 2, and S2 the index variant was classified as a new association if it had a conditional p ≤ 1x10−7 or the P value difference between conditional and unconditional analysis increased by more than two orders of magnitude. Index variants were classified as known if they have previously been reported or the association signal disappeared after conditioning on the variant of a previously reported locus.

103PRS were created for all osteoarthritis phenotypes with 3 different approaches by using 2 different software, PRSice2 which has a P value thresholding shrinkage strategy and LDpred that uses a Bayesian approach to polygenic risk scoring. First, we recreated summary statistics from the main sex-combined meta-analyses excluding arcOGEN samples (Table S1), in order to use the largest possible discovery sample for calculating the weights. PRS were created for the arcOGEN cohort individuals for all available phenotypes (Table S4) using raw genotype data in the software PRSice2 with the binary trait settings. The P value thresholds ranging from 1.3x10−8 to 1.0x10−4 with LD clumping parameters of R2 = 0.1 over 1Mb windows and 10,000 permutations to account for the overfitting. As arcOGEN doesn’t have data for spine osteoarthritis, hand osteoarthritis, thumb osteoarthritis and finger osteoarthritis, we calculated PRS in the UK Biobank samples. For the weights we used the main sex combined meta-analysis excluding UK Biobank individuals. We performed the analysis in PRSice2 using P value thresholds from 1.3x10−8 to 1.0x10−4 with LD clumping parameters of R2 = 0.1 over 1Mb windows and the binary trait setting. In our third approach, deCODE was used as the target dataset since all 11 phenotypes are available and the loss of power in the base dataset wouldn’t be as great as if we had excluded the UKBB dataset. We used PRS analyses for one osteoarthritis trait to investigate its predictive power for another osteoarthritis trait (Table S4). We used effect estimates based on meta-analysis for osteoarthritis excluding deCode. The risk scores were calculated using genotypes for about 600,000 well-imputed autosomal markers. We estimated LD between markers using 4,000 phased Icelandic samples and used this LD information to calculate adjusted effect estimates using LDpred. We created several PRS assuming different fractions of causal markers (the P parameter in LDpred), and selected the PRSs that best predicted the trait itself. The correlation between the PRS and traits was calculated using logistic regression in R (v3.5) adjusting for principal components, sex and year of birth by including them as covariates in the analysis. Variance explained is estimated using Nagelkerke R2. We binned the UK Biobank individuals based on their PRS into deciles and we calculated the odds ratio (OR), 95% CI and P value (Fisher’s Exact test) for individuals in the top decile compared to the bottom decile (Table S4).

104We generated PRS in the genome-wide meta-analysis excluding the deCODE dataset and used univariate linear regression to test the predictiveness of the scores on age at joint replacement in the deCODE cohort (Table S4).

105To investigate PRS predictability of patient strata, we assigned the high confidence effector genes into 6 broad areas of osteoarthritis biological action; skeletal development, joint degeneration, immune function and inflammation, neuronal function and development, muscle function and adipogenesis (Figure 3A; Table S13). For each group the lead SNV for each effector gene member was used to construct a PRS. A meta-analysis without UK Biobank included was used as the base data and UK Biobank was used as the target using PRSice2 (Table S4).

106We investigated associations between osteoarthritis PRS and bone mineral density (BMD) and body mass index (BMI) trajectories in Avon Longitudinal Study of Parents and Children (ALSPAC) study, by analyzing 24,844 BMI, total body less head BMD (TBLH BMD), total body fat mass (TBFM) and total body lean mass (TBLM) observations from 6,263 individuals.

107BMI, total body less head BMD (TBLH BMD), total body fat mass (TBFM) and total body lean mass (TBLM) were measured at the following clinics for adolescents: a) Age 9 (mean age 9.9), b) Age 11 (mean age 11.8), c) Teen focus 2 (mean age 13.8), d) Teen focus 3 (mean age 15.5), e) Teen focus 4 (mean age 17.8) and f) Focus at 24 (mean age 24.5). TBFM and TBLM indices were generated by dividing TBFM or TBLM by height2. Mean BMI, BMD, TBFM and TBLM at each clinic are presented in Table S4. SNVs associated with osteoarthritis at genome-wide significance were included in the scores. Unweighted polygenic risk scores were generated by summing the dosage of the osteoarthritis risk-increasing alleles. Weighted polygenic risk scores were generated by multiplying the dosage of the risk-increasing alleles by the log odds and summing across all alleles. PRS were generated for the SNPs associated with osteoarthritis at any site, hip osteoarthritis and knee osteoarthritis.

108Associations between PRS and BMI trajectories were determined using mixed-effects linear spline regression modeling. 5 knot points were generated at the mean age at each clinic starting at mean age 11.8. PRS-by-spline interaction terms were included to determine if PRS affects the rate of BMI change between time points. Models were adjusted for sex and 4 principal components (PCs). Age (centered at 9.9, the mean age at the first clinic) and individual ID were included as random effects. Analyses were repeated stratified by gender.

109As a sensitivity analysis, analyses were repeated with DXA scans coded as having artifacts excluded, as well as DXA scans coded as having major positioning errors.

110We analyzed six different exposures (weighted and unweighted PRS for osteoarthritis at any site, hip osteoarthritis and knee osteoarthritis) and five different outcomes (BMI, TBLH BMD, TBFMi, TBLMi and height). However, the weighted and unweighted PRS are highly correlated and BMD/TBFMi/TBLMi and height are all components of BMI. Therefore, our corrected P value threshold was determined as 0.05/3 = 0.017. Removing observations extracted from DXA scan images with artifacts or positioning errors did not affect conclusions. Diagnostic tests did provide some evidence of heteroskedasticity and non-normality of residuals for BMI and TBFMi, but log-transforming these outcomes did not alter results or improve residuals. The study website contains details of all the data that is available through a fully searchable data dictionary and variable search tool (http://www.bristol.ac.uk/alspac/researchers/our-data/).

111Results from all independent lead SNVs (n = 100) across all osteoarthritis phenotypes were extracted from the full meta-analysis results. All OR were calculated on the minor allele (allele frequency < 50%), and SNVs with a MAF < 1% were excluded (n = 6). For all the remaining SNVs (n = 94) the OR for each osteoarthritis phenotype GWAS was plotted in a heatmap, together with the corresponding association P (Figure 2). All figures were plotted using R and adjusted for publication quality using Adobe illustrator.

112We also created three classification groups: 0 = Weight bearing joints only (hip and/or knee, knee, hip, total joint replacement, total knee replacement, total hip replacement and spine), 1 = Both, 2 = Non-weight bearing joints only (hand, finger, thumb). Osteoarthritis at any site wasn’t included in this analysis as it wasn’t clear which osteoarthritis subphenotype was leading the signal. In Table 2 and Figure 2 each of the 100 independent genome wide significant SNVs was assigned to the above groups only if it wasn’t nominally significant (p > 0.05) for any of the other phenotypes in the other classification groups, resulting in 86 SNVs to be further analyzed.

113We estimated the genetic correlation between osteoarthritis traits and secondary traits using the cross-trait LD Score regression method as implemented in LDHub (https://github.com/bulik/ldsc). We used results for about 1.1 million well imputed variants, and for LD information we used pre-computed LD scores for European populations (https://alkesgroup.broadinstitute.org/LDSCORE/eur_w_ld_chr.tar.bz2). LD scores for the East Asian populations could not be calculated as the LD Score method requires a sample size of > 4000 individuals.

114To avoid bias due to overlapping samples, we calculated the genetic correlation between equally sized, non-overlapping subgroups of the sample sets from the meta-analysis of Icelandic, Norwegian and USA samples and UK, Dutch, Estonian and Greek samples (Table S6). The results of the two analyses were subsequently meta-analyzed. For genetic correlations with other traits, we calculated the genetic correlation between a meta-analysis of UK, Dutch, Estonian and Greek samples and the Icelandic GWAS summary statistics for each secondary trait, and also between a meta-analysis of Icelandic, Norwegian and USA samples and UKBB GWAS summary statistics for each secondary trait. The results of the two analyses were subsequently meta-analyzed. We also analyzed the genetic correlation between the osteoarthritis subtypes, we split the sample-sets of the meta-analysis in two equally sized groups and performed LD score regression between the two groups for each subtype (Table S6).

115For each independent signal we included all variants within 1Mb of the index variant. In situations where there was more than 1 osteoarthritis signal in the region we used the conditional summary statistics of the meta-analysis conditioned on all other signals. We calculated Wakefield’s asymptotic Bayes factors and we determined the posterior probability of each variant being causal. To produce a 95% credible set of variants we ranked according to posterior probability and included those variants with the highest probability of being causal until the shared probability was at least 95%. Some regions were large therefore we considered only variants in the 95% credible with a posterior probability of causality > 3% (Table S8).

116For cis-eQTL colocalization we used summary statistics of SNPs from 48 human tissues from the GTEx v7 (https://www.gtexportal.org/home/). For each signal and each tissue we included genes that contained at least 1 eQTL (using a threshold of < 5% false discovery rate) in GTEx and that overlapped 100kb either side of our signal. For the colocalization analysis we included all variants in common between the meta-analysis conducted here and the GTEx cis-eQTL analysis with the exception of indels. We used the Bayesian statistical methodology (https://github.com/tobyjohnson/gtx/blob/526120435bb3e29c39fc71604eee03a371ec3753/R/coloc.R) which implements the method of. This method evaluates whether the GWAS and molecular QTL associations best fit a model in which the associations are due to a single shared variant, summarized by the posterior probability (PP). Evidence for colocalization was assessed using the PP4 indicating that there is an association for both traits and they are driven by the same causal variant. A PP4 > 0.8 was considered evidence for colocalization (Table S8).

117Most complex disease risk variants are thought to exert their risk by affecting regulation of expression of a target gene in a tissue and cell-specific context. Osteoarthritis affects multiple tissues within the joint, most notably the cartilage and bone, but there is also evidence for involvement of the synovium, and possibly the muscles and tendons of the joint. To identify possible further osteoarthritis target tissues, we selected all independent genome-wide significant SNVs across osteoarthritis GWAS (n = 100). For each signal we investigated the lead SNV and all fine-mapped SNVs (95% posterior probability) (n = 542) to see if they colocalized with active enhancer histone marks as defined by the ROADMAP epigenomics project (http://www.roadmapepigenomics.org/). All tissue and cell types available in the ROADMAP epigenomics and ENCODE project were used (n = 127). For each lead SNV and their fine-mapped SNVs, the percentage of SNVs located in active enhancer marks was calculated. For the enrichment analysis a background value was calculated (1000 permutations) using 100 random SNVs selected from 1000 Genomes Project (MAF > 0.05). For these 100 random SNVs and all SNV in high LD (R2 = 0.8, LD based on 1000 Genomes Project) their percentage of colocalization with enhancer histone marks was calculated. Once all background permutations were done, an average of all results was taken as the final background values. Enrichment was calculated for each osteoarthritis phenotype and investigated cell type separately, by using the two-proportions Z-test. The significance for enrichment was set to genome-wide significance (p < 1.3x10−8) (Figure S1). As the analysis is highly dependent on the number of variants included (power), significance was only based on enrichment analysis including all 100 independent SNVs across osteoarthritis phenotypes. For the eQTL colocalization using GTEx tissues, we considered the following GTEx tissues as possible osteoarthritis target tissues, based on the tissue specificity analysis: Adipose, Brain, Heart, Lung, Muscle, Nerve, Ovary, Placenta, Skin, Stomach, Cultured fibroblasts, Adrenal gland, and Breast tissue. In the section ‘Amassing evidence to identify effector genes’, genes that had an eQTL colocalization in one of these tissues received an additional scoring point in the lines of evidence (Tables 3 and S10).

118Two-sample Mendelian randomization (MR) was applied to understand the association between plasma proteins on osteoarthritis. In the MR analysis, 1640 proteins in up to 6,000 individuals were treated as the exposure and the 11 osteoarthritis phenotypes as the outcomes. The genetic instruments of the plasma proteins were obtained from, where the conditional independent pQTLs were pooled from 5 recent GWAS of plasma proteins. The genetic instruments were further split into two groups: 1) cis-acting pQTLs within a 500Kb window from each side of the leading pQTL of the protein were used for the MR analysis; 2) trans-acting pQTLs outside the 500kb window of the leading pQTL were designated as trans instruments. For the MR analysis, the meta-analysis summary statistics of osteoarthritis including UK Biobank participants were used as outcomes. We selected a P value threshold of 0.05, corrected for 11 osteoarthritis phenotypes and the number of independent tests, as our threshold for prioritising MR results for follow up colocalization analyses (number of tests = 18.030; p < 2.77x10−6) (Table S7).

119For 28 protein-osteoarthritis associations that survived the multiple testing threshold in the MR analysis, we further conducted colocalization analysis to distinguish causal effects from genomic confounding due to linkage disequilibrium. A colocalization probability more than 80% in this analysis would suggest that the two association signals are likely to colocalize within the test region. Colocalization analysis was applied to both cis and trans pQTLs. For protein and phenotype GWAS lacking sufficient SNP coverage or missing key information (e.g., allele frequency or effect size) in the test region, we conducted a LD check for the sentinel variant for each pQTL against the 30 strongest SNPs in the region associated with the phenotype as an approximate colocalization analysis. R2 of 0.8 between the sentinel pQTL variant and any of the 30 strongest SNPs associated with the phenotype was used as evidence for approximate colocalization (Table S7).

120Nine protein-osteoarthritis associations showed reliable MR and colocalization evidence (Bonferroni corrected, p < 2.77x10−6 and colocalization probability > 70%) for a total of six proteins on seven osteoarthritis phenotypes (Table S7). Of the eight protein quantitative trait loci (pQTLs) used as genetic predictors of these six proteins, five were in strong LD with missense variants (R2 > 0.8). As missense variants may cause epitope-binding artifacts, we also evaluated the effect of these eight pQTLs on other molecular traits: DNA methylation (meQTL) and gene expression (whole blood eQTLs from eQTLGen and all tissues eQTLs from GTEx). Six of the eight pQTLs were also cis meQTLs and cis eQTLs in the same region, in which four pQTLs are in LD (R2 > 0.3) with the top meQTL and eQTL in the region (Table S7).

121We performed an investigation for all genes within 1Mb (upstream and downstream) of the 100 SNVs in preserved and lesioned cartilage and subchondral bone samples from the same donor were obtained from the Research in Articular osteoArthritis Cartilage (RAAK) study consisting of patients with osteoarthritis who underwent joint replacement surgery due to an end-stage disease. Total RNA from articular cartilage and subchondral bone was isolated using QIAGEN RNeasy Mini Kit (QIAGEN, GmbH, Hilden, Germany). Paired-end 2 × 100 bp RNA-sequencing (Illumina TruSeq RNA Library Prep Kit, Illumina HiSeq2000 and Illumina HiSeq4000) was performed. Strand specific RNA-Seq libraries were generated which yielded a mean of 20 million reads per sample. More details on mapping and quality control (QC) from cartilage are previously described. Methods of subchondral bones RNA sequencing have been previously described. After QC, 35 paired cartilage samples (N = 70) and 24 paired subchondral bone samples (18 paired knee and 6 paired hip samples) remained for further differential expression analysis. Normalization and statistical framework was performed using the DESeq2 v1.20 R package. A general linear model (GLM) assuming a negative binomial distribution was applied followed by a paired Wald-test between preserved and lesioned osteoarthritis cartilage and subchondral bone samples. Benjamini-Hochberg multiple testing corrected P values with significance cut-off of 0.05 are reported as False Discovery Rate (FDR).

122We performed an in-silico investigation for all genes within 1Mb (upstream and downstream) of the 100 SNVs in 38 individuals for which differential expression and differential abundance was available, as described previously. Briefly: DNA, RNA, and protein was extracted from matched intact and degraded cartilage samples from 38 patients undergoing total joint replacement surgery: 29 knee and 9 hip osteoarthritis patients. From each patient, 2 paired cartilage samples were taken; a sample with a low Osteoarthritis Research Society International (OARSI) grade which signifies healthy or macroscopically intact cartilage tissue (intact) and a sample with a high OARSI grade, which denotes highly degenerated (degraded) cartilage tissue. All patients provided full written informed consent before participation. The human biological samples were sourced ethically and their research use was in accord with the terms of the informed consents under an Institutional Review Board (IRB)- or Ethics Committee (EC)-approved protocol. Proteomics analysis was performed on intact and degraded cartilage samples from 24 individuals and gene expression analysis on samples from all 38 patients. For the proteomics we performed LC-MS analysis on the Dionex Ultimate 3000 UHPLC system coupled with the Orbitrap Fusion Tribrid Mass Spectrometer. Abundance values were normalized by the sum of all protein abundances in a given sample, then log2-transformed and quantile normalized. We restricted the analysis to 3917 proteins that were quantified in all samples. Differential abundance was performed using a within-individual paired sample design in limma in R. RNA sequencing was performed on the Illumina HiSeq 2000 (75bp paired-end read length). We determined the gene-level counts from transcript-level quantification using salmon 0.8.219 with GRCh38 cDNA assembly release 87 downloaded from Ensembl. Limma-voom was used to remove heteroscedasticity from the estimated expression data. We tested genes for differential expression using limma in R (with lmFit and eBayes), based on a within-individual paired sample design. For the proteomics and RNA sequencing significance was defined at 5% Benjamini-Hochberg FDR to correct for multiple testing.

123For the cartilage RNaseq differential analysis investigation we combined data from the two different sources (detailed above): 35 paired samples from the RAAK study and 38 paired samples from a UK cohort. We only considered genes that were significantly (FDR < 0.05) differentially expressed in both studies (RAAK cohort and UK cohort) with the same direction of effect (Table S9).

124We performed an in-silico investigation for all genes within 1Mb (upstream and downstream) of the 100 SNVs for gene expression in subchondral bone compared to intact cartilage. Briefly, knee joint samples were collected in 11 Han Chinese patients from the Taiwan OA cohort that had undergone total knee replacement surgery (TKR). The subchondral bone tissues underneath the intact and eroded cartilage were obtained as previously described. RNA was extracted as described. RNA (400ng) per sample were used for cRNA synthesis and amplification. Cyanine 3-labeled cRNA was then purified and hybridized to Agilent whole human genome 44 k microarray chips (Agilent Technologies) according to the manufacturer’s instructions (Agilent Technologies, Santa Clara, CA). The array signal intensities were analyzed by the Agilent GeneSpring GX software (version 11.5; Agilent Technologies). Gene expression values were normalized using quantile normalization; probes with low signal intensities were excluded by setting the filter above 32. The normalized values were log transformed and compared using the t test. Differentially expressed genes were defined at ≥ 2 fold-change with Benjamini-Hochberg corrected p ≤ 0.05 (Table S9).

125For the subchondral bone RNaseq differential analysis investigation we combined data from the two different sources (detailed above): 24 paired samples from the RAAK study and 11 Han Chinese patients from the Taiwan OA cohort. We only considered genes that were significantly (FDR < 0.05) differentially expressed in both studies (RAAK cohort and UK cohort) with the same direction of effect (Table S9).

126To cross-reference all independent lead SNVs across all osteoarthritis phenotypes with many other traits and diseases. We queried the PhenoScannerV2 (http://www.phenoscanner.medschl.cam.ac.uk/) database through the web-based tool on 15/02/2021 and at the time it contained > 5,000 genotype-phenotype association datasets from large-scale genetic association studies. Only associations with P value less than the genome wide significance threshold (p < 5x10−8) were included in the analysis and proxies were not requested.

127We examined the druggability status for the 637 prioritized genes (Table S10), using the druggable gene set as defined by. The approach to define the Druggable Genome extended the concept of Mendelian randomization studies for drug development from individual targets to the whole genome. This druggable genome contained 4,479 genes and it was divided in three tiers of druggable gene sets corresponding to position in the drug development pipeline. Tier 1 (1427 genes) included efficacy targets of approved small molecules and biotherapeutic drugs as well as clinical-phase drug candidates. This tier incorporated the targets of approved drugs (licensed drugs) and drugs in clinical development. Proteins that are targets of approved small-molecule and biotherapeutic drugs were identified using manually curated efficacy target information from release 17 of the ChEMBL database. Tier 2 was composed of 682 genes encoding targets with known bioactive drug-like small-molecule binding partners as well as those with ≥ 50% identity (over ≥ 75% of the sequence) with approved drug targets. This tier incorporated proteins closely related to drug targets or with associated drug-like compounds. Proteins closely related to targets of approved drugs were identified through a BLAST search (blastp) of Ensembl peptide sequences against the set of approved drug efficacy targets identified from ChEMBL previously. Tier 3 contained 2370 genes encoding secreted or extracellular proteins, proteins with more distant similarity to approved drug targets, and members of key druggable gene families not already included in tier 1 or 2 [G protein (heterotrimeric guanine nucleotide–binding protein)–coupled receptors (GPCRs), nuclear hormone receptors, ion channels, kinases, and phosphodiesterases].This tier was further subdivided to prioritize those genes that were in proximity (±50 kbp) to a GWAS SNP from GWAS catalog and had an extracellular location (Tier 3A). The remainder of the genes was assigned to Tier 3B. To test if there is enrichment of genes with supporting evidence in the druggable genome database, we conducted a binomial test (N = 205 prioritized genes in the database, N = 637 prioritized genes in total, N = 4479 genes included in the database out of 20330 coding genes tested).

128To glean further insight into the detailed structured information about drugs and drug targets of the identified effector genes in tier 1 or 2 of the druggable genome database, we integrated information from the DrugBank online database (https://www.drugbank.com), accessed between 3–10 March 2021 (version 5.1.8, released 2021-01-03). Of the 71 tier 1 genes, 58 (23 newly associated with osteoarthritis) have info in DrugBank, of which 18 are high-confidence genes (7 newly associated with osteoarthritis).

129We assessed pathway and gene set signal enrichment for the 637 genes with at least one line of evidence and, separately, for the 77 high-confidence putative effector genes (Table S13). Pathway and gene set enrichment analysis were performed using the ConsensusPathDB-human enrichment software (http://cpdb.molgen.mpg.de/). We used the Wikipathways enrichment and the Gene Ontology Biological processes gene set. Enriched gene sets required a minimum of 5 genes to overlap with the examined gene set. The significance threshold was set at FDR < 0.05.

130We also performed 3 additional stratified pathway analysis on (I) all SNVs associated with only weight bearing joints (165 genes), (II) all SNVs associated with only non-weight bearing joints (24 genes), and (III) all SNVs associated with both weight bearing and non-weight bearing joints (155 genes) by using the Gene2Func function of FUMA (Functional Mapping and Annotation of Genome-Wide Association Studies, https://fuma.ctglab.nl/). Analysis was performed as described in using the following settings: all known genes and transcripts were included in the background gene set, we included the MHC-region, a minimum of 5 genes needed to overlap with the examined gene sets and the significance threshold was set at FDR < 0.05 (Table S13).

131The members of arcOGEN Consortium are John Loughlin, Nigel Arden, Fraser Birrell, Andrew Carr, Panos Deloukas, Michael Doherty, Andrew W. McCaskie, William E.R. Ollier, Ashok Rai, Stuart H. Ralston, Tim D. Spector, Ana M. Valdes, Gillian A. Wallis, J. Mark Wilkinson, and Eleftheria Zeggini.

132The members of the HUNT All-In Pain Consortium are Amy E. Martinsen, Cristen Willer, Egil Andreas Fors, Ingunn Mundal, Knut Hagen, Kristian Bernhard Nilsen, Marie Udnesseter Lie, Sigrid Børte, Ben Brumpton, Jonas Bille Nielsen, Lars G. Fritsche, Wei Zhou, Ingrid Heuch, and Kjersti Storheim.

133The members of the ARGO-Athens Consortium are Eleni Zengini, George Alexiadis, Evangelos Tyrpenou, Athanasios Koukakis, Dimitrios Chytas, Dimitrios Stergios Evangelopoulos, Chronopoulos Efstathios, Spiros Pneumaticos, Vasileios S. Nikolaou, J. Mark Wilkinson, George C. Babis, Konstantinos Hatzikotoulas, and Eleftheria Zeggini.

134The members of the ARGO-Larissa Consortium are Konstantinos Malizos, Lydia Anastasopoulou, Aspasia Tsezou, Eleni Zengini, J. Mark Wilkinson, Konstantinos Hatzikotoulas, and Eleftheria Zeggini.

135The members of the Regeneron Genetics Center Consortium are Goncalo Abecasis, Aris Baras, Michael Cantor, Giovanni Coppola, Andrew Deubler, Aris Economides, Luca A. Lotta, John D. Overton, Jeffrey G. Reid, Alan Shuldiner, Katia Karalis, Katherine Siminovitch, Christina Beechert, Caitlin Forsythe, Erin D. Fuller, Zhenhua Gu, Michael Lattari, Alexander Lopez, Thomas D. Schleicher, Maria Sotiropoulos Padilla, Louis Widom, Sarah E. Wolf, Manasi Pradhan, Kia Manoochehri, Xiaodong Bai, Suganthi Balasubramanian, Boris Boutkov, Gisu Eom, Lukas Habegger, Alicia Hawes, Olga Krasheninina, Rouel Lanche, Adam J. Mansfield, Evan K. Maxwell, Mona Nafde, Sean O’Keeffe, Max Orelus, Razvan Panea, Tommy Polanco, Ayesha Rasool, William Salerno, Jeffrey C. Staples, Dadong Li, Deepika Sharma, Ilanjana Banerjee, Jonas Bovijn, Adam Locke, Niek Verweij, Mary Haas, George Hindy, Tanima De, Parsa Akbari, Olukayode Sosina, Manuel A. R. Ferreira, Marcus B. Jones, Jason Mighty, Michelle G. LeBlanc, and Lyndon J. Mitnaul.