Prioritization of PLEC and GRINA as Osteoarthritis Risk Genes Through the Identification and Characterization of Novel Methylation Quantitative Trait Loci

Abstract

Objective: To identify methylation quantitative trait loci (mQTLs) correlating with osteoarthritis (OA) risk alleles and to undertake mechanistic characterization as a means of target gene prioritization.

Methods: We used genome-wide genotyping and cartilage DNA methylation array data in a discovery screen of novel OA risk loci. This was followed by methylation, gene expression analysis, and genotyping studies in additional cartilage samples, accompanied by in silico analyses.

Results: We identified 4 novel OA mQTLs. The most significant mQTL contained 9 CpG sites where methylation correlated with OA risk genotype, with 5 of the CpG sites having P values <1 × 10^-10 . The 9 CpG sites reside in an interval of only 7.7 kb within the PLEC gene and form 2 distinct clusters. We were able to prioritize PLEC and the adjacent gene GRINA as independent targets of the OA risk. We identified PLEC and GRINA expression QTLs operating in cartilage, as well as methylation-expression QTLs operating on the 2 genes. GRINA and PLEC also demonstrated differential expression between OA hip and non-OA hip cartilage.

Conclusion: PLEC encodes plectin, a cytoskeletal protein that maintains tissue integrity by regulating intracellular signaling in response to mechanical stimuli. GRINA encodes the ionotropic glutamate receptor TMBIM3 (transmembrane BAX inhibitor 1 motif-containing protein family member 3), which regulates cell survival. Based on our results, we hypothesize that in a joint predisposed to OA, expression of these genes alters in order to combat aberrant biomechanics, and that this is epigenetically regulated. However, carriage of the OA risk-conferring allele at this locus hinders this response and contributes to disease development.

Figure 1

Correlation of genotype at rs11780978 with the methylation status of 9 CpG sites within the gene body of PLEC. A, The association between rs11780978 and methylation levels of CpG probes that are present within the region. The x‐axis represents the genomic position of the CpG probes, and the y‐axis represents the Benjamini‐Hochberg–corrected −log₁₀ P value of the correlation between rs11780978 genotype and M value at each CpG probe. Each circle represents a single CpG probe, with the 9 significant associations highlighted in red. The broken line indicates the location of rs11780978. The genes within the region analyzed are indicated below the association plot, with the gene direction indicated by arrows. B, An enlarged image of PLEC, highlighting the location of the 9 significant CpG associations (red) falling into 2 distinct clusters. Each circle represents a single CpG probe. C, The association between genotype at rs11780978 and methylation levels at the 9 significant CpG probes for all 87 samples. The level of methylation at the CpG probes is shown as the β value. Symbols represent individual patient samples. Horizontal lines and error bars show the mean ± SEM.

Figure 2

Association between rs11780978 genotype and methylation for PLEC CpG sites cg19405177 and cg14598846 in the new osteoarthritis patients (n = 104). In cluster 1, the pyrosequencing assay targeting cg19405177 (CpG5) captured 7 additional CpG sites (CpG1–CpG4 and CpG6–CpG8) located between 67 bp upstream and 45 bp downstream of cg19405177. In cluster 2, the pyrosequencing assay targeting cg14598846 (CpG1) captured 3 additional CpG sites (CpG2–CpG4) located up to 22 bp downstream of cg14598846. P values were calculated using the Kruskal‐Wallis test. The square of the correlation coefficient (r²) values were calculated using a model of standard least squares. Horizontal lines and error bars show the mean ± SEM. n = the number of patients providing data per CpG.

Figure 3

Allelic expression imbalance (AEI) analysis in cartilage from new osteoarthritis patients (n = 104). AEI analysis was conducted for PLEC using transcript single‐nucleotide polymorphism (SNP) rs11783799 (A), for PARP10 using transcript SNP rs11136345 (B), and for GRINA using transcript SNP rs9100 (C). The risk/nonrisk allelic ratio is plotted, with a ratio of <1 indicating decreased expression of the risk allele. Three technical repeats were performed for each patient's DNA sample (black) and cDNA sample (gray). Numbers on the x‐axis refer to the anonymized identification number assigned to patients at recruitment. The mean values for DNA and cDNA are combined (left panels), with results represented as box plots, in which the lines within the box represent the median, the box represents the 25th to 75th percentile, and the whiskers represent the maximum and minimum values (right). P values were calculated using Wilcoxon's matched pairs signed rank test.

Figure 4

Methylation expression quantitative trait locus analysis of PLEC (A) and GRINA (B). A,PLEC log₂ allelic expression imbalance (AEI) ratios and B,GRINA log₂ AEI ratios were plotted against methylation at cg19405177 and its 7 additional CpG sites (cluster 1) and at cg14598846 and its 3 additional CpG sites (cluster 2).