Evolutionary Selection and Constraint on Human Knee Chondrocyte Regulation Impacts Osteoarthritis Risk

Abstract

During human evolution, the knee adapted to the biomechanical demands of bipedalism by altering chondrocyte developmental programs. This adaptive process was likely not without deleterious consequences to health. Today, osteoarthritis occurs in 250 million people, with risk variants enriched in non-coding sequences near chondrocyte genes, loci that likely became optimized during knee evolution. We explore this relationship by epigenetically profiling joint chondrocytes, revealing ancient selection and recent constraint and drift on knee regulatory elements, which also overlap osteoarthritis variants that contribute to disease heritability by tending to modify constrained functional sequence. We propose a model whereby genetic violations to regulatory constraint, tolerated during knee development, lead to adult pathology. In support, we discover a causal enhancer variant (rs6060369) present in billions of people at a risk locus (GDF5-UQCC1), showing how it impacts mouse knee-shape and osteoarthritis. Overall, our methods link an evolutionarily novel aspect of human anatomy to its pathogenesis.

Keywords: ATAC-seq; GDF5; chondrocyte; gene regulation; genetic drift; human evolution; knee; mouse model; natural selection; osteoarthritis.

Figure 1.. Sequence Features of Chondrocyte Epigenetic Profiles

(A) Diagram of ATAC-seq tissues. Numbers indicate replicate-consolidated peak calls before filtering with brain data. Inset: genomic distribution of lifted-over DF (top) and PT (bottom) elements. (B) Per-bp conservation scores for 20 primates (phyloP20ways) aggregated for ATAC-seq sets and compared via Wilcoxon rank-sum test. (C) Overlaps of human acceleration regions and ATAC-seq sets. KC, knee-common; GM, GM12878; DF, distal femur-specific; PT, proximal tibia-specific; EB, embryonic brain. Overlaps shown relative to set size (per bp of sequence) for background (gray) and target (colored) sets. (D) Relationship between human-chimp sequence identity (red line, right y axis) and phyloP20ways score averaged over an ATAC-seq peak (blue line, left y axis) for DF-specific elements. Dashed curve and horizontal lines indicate average phyloP20ways scores for random background region set. Shaded regions reflect low/high %ID regions used in GREAT. Significance codes: not significant (ns), < 0.05 (*), < 0.01 (**), < 1e–5 (***). See also Figure S1 and Tables S1, S2, and S4.

Figure 2.. Sequence Modification of Putative Regulatory Elements

(A) Histograms showing depletion and enrichment of human-chimp base pair alterations intersecting predicted FOXP1 and KLF5 TF motifs, respectively, in DF-specific elements relative to randomized sets (p < 0.01). Counts of altered bases in background (gray) and element (red) sets shown as a fraction of total set size. See Table S3. (B) KLF5 motif example in a DF-specific element intronic to PBX1, a TF regulating HOX expression and endochondral ossification (Capellini et al., 2006; Selleri et al., 2001). Predicted disruption by a human-derived T/G nucleotide change shown. A region of acceleration in a knee element identified upstream of PBX1 (Table S3). Red shading indicates altered base relative to motif logo, blue indicates genomic position of sequence. ATAC-seq regions and phyloPways conservation tracks, UCSC sequence alignment shown. (C) Example of a human-chimp alteration predicted to improve a CTCF motif within a PT-specific element intronic to NCOR2, involved in skeletal biology (Blake et al., 2017). Format similar to (B).

Figure 3.. Intra-species Variation in Regulatory Elements and Human Variation in Knee Morphology

(A) Counts of common human variants per bp of sequence for element sets compared to random region sets along with other genomic features; labels correspond to Table S5. (B) Common human variants intersecting elements in the elbow (left) and knee (right) specific sets were counted and compared across sets. (C) Comparisons of chimp sequence diversity across knee-specific sets. (D–F) Common variants in human, chimp, and gorilla intersecting a given ATAC-seq peak counted for all variable sequences in knee-specific sets, expressed as “SNPs per sequence”—mean/median values shown in dashed/bold lines, respectively. (D) Distal femur-specific (DF). (E) Proximal tibia-specific (PT). (F) Knee-common-specific (KC). (G) Measurements of human medial/lateral tibia and femur via MRI dataset. Volume of cartilage (left) and total area of subchondral bone (right) for medial/lateral portions of both bones were measured across all subjects. Individual points are plotted alongside an outlined density curve, quartiles indicated in dashed lines. Significance codes: not significant (ns), < 0.05 (*), < 0.01 (**), < 1e–5 (***). See also Figure S2 and Tables S5 and S6.

Figure 4.. Human Variation in Knee Elements and Its Impacts on OA Risk

(A) Enrichments for OA risk variants in general knee, knee-specific, human E59 hind limb, and brain region sets, along with human E47 and BMDC H3K27ac ChIP-seq regions. Calculated Z score enrichments over randomized sets shown as −log (p value); red line indicates significance threshold (p < 0.05). (B) Distribution of predicted motifs intersected by OA risk variants, counted per TF. Significantly enriched factors are indicated—FOXP1 significant following p value correction (p < 0.05). (C) Distribution of predicted motifs intersecting OA risk variants in knee elements for a set of chondrogenesis-related TFs. (D) Overlaps of region sets and signals of recent selection calculated for general knee, knee-specific, brain, and human H3K27ac ChIP-seq region sets, along with OA risk-variants and the regions intersecting them. Clustered set overlaps also shown. Hypergeometric tests represented as −log (p value), with sign denoting enrichment/depletion; red lines indicate respective significance cut-offs (p < 0.05). (E) Top: formalized model for the role of evolutionary history in modern heritable OA risk. Ancient selection acting on ancestral sequence diversity in regulatory elements establishes a derived knee configuration, which is subsequently maintained through ancient purifying selection (i.e., functional constraint). More recently, genetic drift, in combination with antagonistic selection for other traits, increases the frequency of alternative alleles in functionally constrained sites. Bottom: the presence of moderate mutational load in highly constrained elements (e.g., prox. tibia elements), or high mutational load in less-constrained elements (e.g., dist. femur elements), stand to disrupt knee homeostasis and promote pathology risk, while low mutational loads (or low sequence constraint) are more tolerated (i.e., harbor lesser risk of pathology). (F) The number of alternative alleles falling in chondrocyte regulatory elements counted per-individual for the 1KG3 population (blue) and the OAI patient cohort (red), shown as density distributions with mean values (dashed lines); significance bar indicates Student’s t test result (p < 0.05). See also Figure S3 and Tables S3, S7, and S8.

Figure 5.. Functional Characterization of the GDF5 Locus and R4 Enhancer in the Mouse

(A) UCSC Genome Browser view of human GDF5 locus with intersections of OA variants, knee ATAC-seq regions (for human and mouse tissue), GDF5 enhancers, and H3K27ac ChIP-seq signals from human embryonic limbs (E47) and BMDCs. Rs4911178 (left red line) and rs6060369 (right red line) overlap functional knee sequences. (B) R4-driven lacZ expression in inferior-most distal femur (DF) and superior-most proximal tibia (PT) (left) tissues and in adult distal femora (right two images). AC, articular cartilage; FN, femoral notch; IC, inferior condyles; CL, cruciate ligaments; TG, trochlear groove. Scale bars, E14.5 = 250 μm; adult = 1 mm. (C) P30 knee anatomy in C57BL/6J R4 null mice. #, same trends for medial condyle. R4^+/+ (wild-type [WT]), R4^−/− homozygous (HOM). (D) 1-year knee anatomy. (E) OARSI scores on WT and HOM knees at P30 and 1 year. Triangles, heterotopic ossification. (F) 3D renditions (top) and histology (bottom) of WT/HOM knees with minimum (blue), mean (green), and maximum (red) OARSI scores. Heterotopic ossification observed in HOM knees with highest OARSI scores and most cartilage damage. Scale bars, 50 μm. (G) X-ray 3D cartilage scanning of WT/HOM distal femur condyles (top) and proximal tibia platforms (bottom) at 1 year, showing OA lesions (white arrows). Scale bars, 1 mm. (H) Linear regressions between knee morphology and OARSI score at 1 year. In (C)–(E) colored dots correspond to specimens shown in (F). Significance code: <0.05 (*). See also Figures S4 and S5 and Table S9.

Figure 6.. Functional Characterization of the R4 Enhancer in Human and Mouse Chondrocytes

(A) Expression by qRT-PCR of CEP250 (not significant [n.s.]), GDF5 (p < 0.005), and UQCC1 (n.s.) in human T/C-28a2 chondrocytes lacking R4. (B) Expression of CEP250 (n.s.), GDF5 (p < 0.05), and UQCC1 (n.s.) in cells lacking 41 bp containing rs6060369 in R4. (C) In vitro reporter analyses of R4-driven luciferase activity comparing constructs with the OA risk “T” or non-risk “C” variant in T/C-28a2 cells. (D) PITX1 ChIP on a R4 sub-region containing rs6060369 in T/C-28a2 cells showing input (left) and pull-down (right) in image. (E) UCSC Genome Browser view of mouse R4 corresponding to knee ATAC-seq region and PITX1 binding via ChIP-seq (Infante et al., 2013).

Figure 7.. Functional Characterization of the rs6060369 OA Risk Allele in Humanized Replacement Mice

(A) ASE assays on E14.5 distal femur chondrocytes from C57BL/6J/129SVJ R4^{rs6060369-T, rs6060369-A} replacement mice (right, p = 0.0005) and C57BL/6J/129SVJ R4 heterozygous (R4^+/−) mice (left, p = 0.001). (B) 3D comparative analysis indicating the locations of largest morphological differences between (left) WT R4^rs6060369-A, rs6060369-A and HOM R4^{rs6060369-T, rs6060369-T} hind limbs, as well as between (right) R4^+/+ and R4^−/− hind limbs (zoom-in images focus on inferior distal femur [top], superior proximal tibia [bottom]). Areas with largest variations are highlighted in red (WT > HOM) and dark blue (WT < HOM), with minimal variation in green/yellow. (C) mCT measurements of indicated features in base pair replacement and R4 mice at postnatal days P30, P56, and 1 year. See also Figure S4 and Table S9.