Genetics of osteoarthritis

Abstract

Osteoarthritis genetics has been transformed in the past decade through the application of large-scale genome-wide association scans. So far, over 100 polymorphic DNA variants have been associated with this common and complex disease. These genetic risk variants account for over 20% of osteoarthritis heritability and the vast majority map to non-protein coding regions of the genome where they are presumed to act by regulating the expression of target genes. Statistical fine mapping, in silico analyses of genomics data, and laboratory-based functional studies have enabled the identification of some of these targets, which encode proteins with diverse roles, including extracellular signaling molecules, intracellular enzymes, transcription factors, and cytoskeletal proteins. A large number of the risk variants correlate with epigenetic factors, in particular cartilage DNA methylation changes in cis, implying that epigenetics may be a conduit through which genetic effects on gene expression are mediated. Some of the variants also appear to have been selected as humans adapted to bipedalism, suggesting that a proportion of osteoarthritis genetic susceptibility results from antagonistic pleiotropy, with risk variants having a positive role in joint formation but a negative role in the long-term health of the joint. Although data from an osteoarthritis genetic study has not yet directly led to a novel treatment, some of the osteoarthritis associated genes code for proteins that have available therapeutics. Genetic investigations are therefore revealing fascinating fundamental insights into osteoarthritis and can expose options for translational intervention.

Keywords: DNA methylation; Epigenetics; Functional analysis; GWAS; Genetics; SNPs.

Fig. 1

The genomic locations of the SNPs reported to be associated with OA. UTR, untranslated region.

Fig. 2

Example of how an OA GWAS risk SNP can lead to mechanistic insights and an associated OA risk gene. (A) Manhattan plot of an OA GWAS. Each dot represents a SNP. The -log10 P-value represents the significance of each SNP being preferentially carried by OA patients vs healthy controls. The dotted line represents the genome wide significance threshold. SNPs rs1, rs2 and rs3 are all OA significant risk SNPs. (B) Cartilage mQTL analysis for rs1. Each circle represents a CpG site available on the methylation array. The -log10 P-value represents the significance that methylation at a CpG is associated with rs1 genotype. The dotted line represents the statistical significance threshold after multiple testing correction. The red circle represents a CpG site (cg1) significantly associated with rs1 genotype (this is an mQTL). This mQTL analysis should be done for each SNP that is found to be associated with OA (rs2 and rs3). This analysis is usually limited to the CpGs within the physical proximity of each association SNP (e.g., 1 megabase (Mb) region). (C) Arrows representing the location and transcriptional direction of each gene within the region. rs1 is located within GENE3. (D) Linkage disequilibrium (LD) in CEU population from the HapMap project. Image correlates with the topologically associated domains (TADs), shown by each pyramid, within the region. Red represents regions in LD, with each pyramid encompassing a region that is inherited as a block. In this example both rs1 and cg1 are within the same TAD as are GENE2 and GENE3. These genes would be prioritized as candidate OA risk genes for further analysis. (E) LD link r² values between rs1 and each SNP within the region. Each bar represents a SNP and the height corresponds to the r² value between 0 and 1. All SNPs in high LD with rs1 are just as likely to be the functional SNP. (F) UCSC (http://genome.ucsc.edu/) track from GTEx showing all eQTLs operating within the region. Each line represents an eQTL in a tissue type for one of the four genes listed. The position of the bar on the x-axis represents which SNP this operates at. All SNPs with an eQTL are prioritized as functional. (G) UCSC track of ChIP-seq data in chondrocytes and osteoblasts. Green represents transcription sites, yellow enhancer sites and red transcription start sites. Enhancers and transcription start sites are prioritized as functional. (H) Genome wide ATAC-seq data for knee articular cartilage available on UCSC browser. (I) Association between genotype at rs1 and GENE2 expression levels using qPCR. Each dot represents cartilage tissue from a different patient. C allele is associated with increased expression. (J) Violin plot showing association between genotype at rs1 and methylation levels at cg1. The width of the blue violin represents the proportion of patients within that range, the bar the mean and the dotted line the quartile range. C allele at rs1 is associated with decreased methylation at cg1. (K) Allelic expression imbalance (AEI) of GENE2. rs1 C allele is preferentially transcribed compared to the T allele. (L) Association between GENE2 expression levels and cg1 methylation levels, highlighting a methylation-expression QTL (meQTL). Each dot represents cartilage tissue from a different patient. Increased GENE2 expression is associated with decreased cg1 methylation levels. (M) Association between GENE2 AEI and cg1 methylation, again highlighting an meQTL. An increased ratio of rs1 C/T allele transcripts is associated with decreased methylation at cg1. (N) In vitro Lucia reporter analysis with DNA insert containing rs1 and cg1. rs1 C and T alleles are compared to each other with cg1 methylated or unmethylated. C allele and decreased methylation act synergistically to increase enhancer activity. (O) In vitro deletion of cg1 region using CRISPR-Cas9 and subsequent gene expression analysis of GENE2 and GENE3. Each dot represents a biological in vitro replicate. GENE3 expression levels are unchanged whilst GENE2 expression levels are reduced. (P) In vitro methylation and demethylation of cg1 using CRISPR-dCas9 DNMT3a/TET1. On the x-axis, controls are listed by the abbreviation C and cases by the lightning bolt. Methylation changes do not influence GENE3 expression levels. Increased methylation of cg1 decreases GENE2 expression and conversely, decreased methylation increases expression levels. This concludes that methylation at cg1 is only driving GENE2 expression. GENE2 is deemed the target of the OA risk marked by SNP rs1.

Fig. 3

A STRING protein–protein interaction network of OA risk gene products. Yellow nodes indicate proteins encoded by OA risk genes COLGALT2, RUNX2, PLEC, MGP, TGFB1 and GDF5. Blue nodes are proteins encoded by unconfirmed causal genes SMAD3 and CDC5L at OA risk loci. Grey nodes are protein interactors. Nodes were functionally annotated with the top three Gene Ontology (GO) processes: extracellular matrix (ECM) organization, green (P = 5.23 × 10⁻¹⁴); tissue development, orange (P = 5.06 × 10⁻¹³); and cellular response to growth factor stimulus, purple (P = 1.63 × 10⁻¹³). Edge thickness (lines between proteins) indicates the confidence of interactions, based upon experimental evidence, co-expression, databases, and text-mining. STRING website; https://string-db.org/.