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

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.


INTRODUCTION
Osteoarthritis (OA) is a common, age-related disease of synovial joints characterized by the thinning of articular cartilage. This thinning ultimately results in focal loss of cartilage and a fullthickness lesion exposing underlying bone (1,2). Cartilage loss prevents the joint from withstanding normal mechanical load and leads to severely impaired joint function. The clinical effect is a painful chronic morbidity and, in those with hip and knee OA, an increased risk of premature death due to secondary cardiovascular events resulting from reduced physical activity (3,4). Regarding treatment options for hip and knee OA, there are no disease-modifying OA drugs, and arthroplasty of hips and knees is a common, but not a risk-free, procedure. As the population ages, the prevalence of OA increases. Novel treatments are therefore urgently required.
The causes of OA are complex and, with a heritability of >40%, genetic susceptibility is a main driver (5). Genome-wide association studies have revealed that the OA genetic component is polygenic, with the association signals typically having modest odds ratios of <1.5 (6).
Of those OA risk-conferring loci that have so far been functionally investigated, the vast majority mediate their effect by modulating the expression of a nearby gene. Clear examples include the single-nucleotide polymorphisms (SNPs) rs143383, rs225014, and rs3204689, which correlate with differential expression of the OA and non-OA risk alleles in cartilage of the genes GDF5, DIO2, and ALDH1A2, respectively (for review, see ref. 6).
DNA methylation at CpG dinucleotides is an epigenetic process used by the cell to regulate gene expression, and it can amplify or attenuate the impact of an expression quantitative trait locus (QTL) (7)(8)(9). Our group has previously shown this in the context of GDF5, rs143383, and CpG sites located close to the SNP (10). Identifying cis-acting CpG sites, where methylation correlates with an association signal, can therefore offer insight into the mechanism by which the genetic risk at that signal operates. Such CpG sites are known as methylation QTLs (mQTLs). Furthermore, if mQTL CpG sites cluster within or close to a particular gene, their mapping prioritizes that gene for further investigation as a plausible target for the association signal (11).
We have shown the potential utility of this approach in our previous analysis of 16 OA loci, in which we correlated cartilage DNA methylation with OA association signals using genome-wide DNA methylation and genotyping array data (12). In this current study, we extended our investigation to 18 novel OA association signals that have been identified in the last 3 years. We found evidence of mQTLs at 4 of these loci, with the most statistically significant effect being on chromosome 8q24.3 at a signal harboring a number of genes. We subjected this signal to in silico and in vitro analyses, which highlighted the plectin gene PLEC and the TMBIM3 (transmembrane BAX inhibitor 1 motif-containing protein family member 3) gene GRINA as the targets of the association. Plectin is a multifunctional cytoskeletal protein that is particularly abundant in tissues subjected to mechanical load and stress, while TMBIM3 is involved in the control of cell death by endoplasmic reticulum (ER) stress. Our study highlights the high frequency of cartilage mQTLs operating on OA risk loci and the capacity for mQTL analysis to prioritize a gene from within a gene-rich region as a target of the association signal. PLEC, GRINA, and their encoded proteins now merit much more detailed analyses in the context of OA genetic risk and cartilage biology.

MATERIALS AND METHODS
Cartilage methylation and genotyping array data. We used cartilage CpG methylation and genotype data that we had previously generated using an Illumina Infinium HumanMethylation450 array and an Illumina HumanOmniExpress array, respectively (12). We had both methylation and genotype data available for a total of 87 patients who had undergone knee or hip joint arthroplasty (57 knee OA patients, 14  OA loci investigation and mQTL analysis. We investigated 18 novel OA risk loci that had been reported as being associated with the disease at a significance level close to or surpassing the genome-wide threshold of P < 5 × 10 −8 (13)(14)(15)(16)(17)(18)(19)(20) ( Table 1). If the OA-associated SNP was directly genotyped on an Illumina HumanOmniExpress array, we utilized those SNP data. If the association SNP was not on the array, we identified and, where possible, used a proxy SNP that was in perfect or high linkage disequilibrium (pairwise r 2 > 0.7) with the association SNP. Proxy SNPs were derived from a candidate list using LDlink's LDproxy tool (https ://analy sisto ols.nci.nih.gov/ LDlin k/) (21) and European population data. Where multiple proxy SNPs were identified, we chose the one with the highest r 2 relative to the association SNP.
For each locus, we covered a 2-Mb region encompassing 1 Mb upstream of and 1 Mb downstream of the association SNP. For each CpG site within the 2-Mb region, linear regression was used to measure the relationship between methylation, in the form of M values, and genotype (0, 1, or 2 copies of the minor allele). For the purpose of this analysis, we defined any genotype-methylation correlations identified within this 2-Mb region to be a cis mQTL. CpG sequences that were directly modified by SNPs were not included in the analysis. We used age, sex, and sample type (OA knee, OA hip, and femoral neck fracture) as covariates. All mQTL calculations were performed using Matrix eQTL software (22), which implements a false discovery rate (FDR) estimation based on the Benjamini-Hochberg FDR procedure (23) and accounts for the number of tests performed. Methylation β values were used for the purpose of visualization.
RNA-sequencing (RNA-Seq) analysis. The expression of genes of interest was assessed using RNA-Seq data that we had previously generated from the cartilage of 10 patients with hip OA and 6 patients with femoral neck fracture (24,25). These patients are unrelated to the 87 patients used in the mQTL analysis. Details regarding the 16  Genotyping of SNPs. SNPs were genotyped by pyrosequencing (rs11783799 and rs11136345), by restriction fragment length polymorphism (RFLP) analysis (rs11780978 and rs7819099), or by a real-time SNP genotyping assay (rs9100) (catalog no. 4351379; ThermoFisher Scientific). The rs11136336 SNP is in perfect linkage disequilibrium (r 2 = 1.0) with rs11783799; the rs11783799 assay was therefore used to genotype rs11136336. Pyrosequencing and RFLP polymerase chain reactions (PCRs) were performed using a G-Storm GS4 Q4 Quad Block Thermal Cycler (Somerton Biotechnology) and the primers listed in Supplementary Table 2 (available on the Arthritis & Rheumatology web site at http://onlin elibr ary.wiley.com/doi/10.1002/art.40849/ abstract). Pyrosequencing assays were designed using PyroMark assay design software 2.0 (Qiagen), and the sequencing was performed using a PyroMark Q24 Advanced platform (Qiagen) with the recommended kit, following the instructions of the manufacturer. The RFLP-digested fragments were separated by electrophoresis through a 3% agarose gel and visualized using ethidium bromide staining. The rs9100 real-time genotyping assay was run on a QuantStudio 3 (Applied Biosystems), using TaqPath ProAmp * MAF = minor allele frequency; Chr. = chromosome. † For single-nucleotide polymorphisms (SNPs) not present on the HumanOmniExpress array, a proxy SNP that had the highest linkage disequilibrium with the association SNP and that was on the array was used to infer genotypes. The threshold was set at r 2 > 0.7. Thirteen of the association SNPs required a proxy and, at the threshold, a proxy was available for 9 of the 13 SNPs (rs2820436, rs2862851, rs10471753, rs3850251, rs496547, rs4764133, rs754106, rs2521349, and rs6516886), but a proxy was not available for 4 of the 13 SNPs (rs2236995, rs11335718, rs833058, and rs864839). ‡ Stratum highlights the joint with which the signal shows association from the original genetic study. § This value shows the number of CpG probes from the Illumina Infinium HumanMethylation450 array that are present within the 2-Mb region surrounding the association SNP.
MasterMix (ThermoFisher) following the instructions of the manufacturer.
Complementary DNA synthesis and AEI analysis. Complementary DNA was synthesized from 1 μg of cartilage total RNA using a SuperScript First-Strand synthesis system (Invitrogen) and random hexamers following the standard protocol of the manufacturer after an initial 30-minute treatment with 1 unit of Turbo DNase (Invitrogen) at 37°C.
AEI at transcript SNPs rs11783799 and rs11136345 was quantified by pyrosequencing, using the methodology described in the above genotyping section and the same primers. The sequences were generated automatically, and an output of allelic ratio was produced using PyroMark Advanced software. AEI at transcript SNP rs9100 was quantified using the above real-time SNP genotyping assay. The use of such assays for AEI has been described previously (30). For each cartilage cDNA and DNA sample, PCRs were formed in triplicate. Samples were excluded from the analysis if the values between the PCR replicates differed by >5% for the pyrosequencing assays, or by >0.5 cycle thresholds (C t ) for the genotyping assay. The respective cDNA and DNA were analyzed concurrently, and allelic expression of cDNA was normalized to its corresponding DNA.
Targeted CpG methylation analysis. Methylation analysis of CpG sites cg19405177 and cg14598846 was performed by pyrosequencing using the methodology described above and the PCR and sequencing primers listed in Supplementary Chromatin interactions. The WashU Epigenome Browser (31) was searched to identify long-range chromatin interactions extending from the CpG sites within the PLEC locus. All publicly available Hi-C and long-range chromatin interaction data sets were loaded for all cell types with available data within the genomic region chromosome 8:145,000,040-145,068,742. The region was searched visually to identify interactions stemming from either of the 2 clusters of CpG sites. The human breast cancer cell line MCF7 and the chronic myeloid leukemia cell line K562 interaction schema represent protein factor-mediated chromatin interaction data measured by RNA polymerase II chromatin interaction analysis using paired-end tag sequencing (ChIA-PET) data (32). These data were produced as part of the ENCODE project (https://www.encodeproject.org/). Data from the human lymphoblastoid cell line GM12878 were produced by RNA CCCTCF-binding factor ChIA-PET (33).

RESULTS
Identification of OA mQTLs. Table 1 provides details of the 18 loci investigated. For 5 of the 18 loci, the SNP that had been used to identify the association signal had been genotyped on the HumanOmniExpress array. For each of the remaining 13 loci, we searched for a proxy SNP (r 2 > 0.7). This was successful for 9 loci, but for the remaining 4 loci there was no proxy (loci 4, 5, 9, and 16). Therefore, these 4 loci were not studied further. In total, we analyzed 4,895 CpG sites across the 14 loci. We assessed correlations in all 87 samples combined for knee OA, hip OA, and femoral neck fracture. This analysis identified 4 SNPs that correlated with methylation (FDR P < 0.05): locus 7, SNP rs10471753 and 1 CpG site; locus 11, SNP rs11780978 and 9 CpG sites; locus 14, SNP rs4764133 and 1 CpG site; and locus 18, SNP rs6516886 and 6 CpG sites ( Table 2 and Supplementary Table 4 Many OA genetic association signals have only been discovered following stratification by joint and/or sex (5,6). We repeated the mQTL analysis of all 14 loci using such stratification, but did not identify any further genotype-methylation correlations. Furthermore, to assess whether the identified mQTLs were more specific to OA than femoral neck fracture, we repeated the mQTL analysis by removing the 16 samples from the patients with femoral neck fracture. This analysis did not identify any OA-specific or femoral neck fracture-specific genotype-methylation correlations. In summary, we identified 4 OA susceptibility loci demonstrating significant mQTLs, and these are discussed in more detail below.

Locus 7.
Genotype at rs10471753 correlated with methylation of 1 CpG site, cg25008444, which is located 233 kb from the association SNP. This CpG site is within the gene body of PIK3R1, coding for phosphatidylinositol 3-kinase regulatory subunit ɑ (Supplementary Figure 1A, Figure 1B).

Locus 11.
Genotype at rs11780978 correlated with methylation of 9 CpG sites located within an interval of only 7.7 kb and positioned between 25.7 kb and 33.4 kb from the association SNP. The 9 CpG sites are all within the gene body of PLEC, which encodes plectin ( Figure 1A). They form 2 clusters as follows: cg19405177, cg20784950, and cg01870834 are located within an interval of 1.4 kb, followed by a gap of 5.3 kb (which contains 6 CpG sites covered by the array and not demonstrating an mQTL) and then 6 CpG sites (cg07427475, cg02331830, cg04255391, cg14598846, cg23299254, and cg10299941), which are located within an interval | 1289 of 1.0 kb ( Figure 1B). For each of the 3 CpG sites within the first cluster, the OA risk-conferring A allele of rs11780978 correlated with higher methylation, while for each of the 6 CpG sites within the second cluster, the A allele correlated with lower methylation ( Figure 1C).

Locus 14.
Genotype at rs4764133 correlated with methylation at 1 CpG site, cg20917083, which is located 50 kb from the association SNP. This CpG is within an intron of ARHGDIB, which codes for Rho GDP dissociation inhibitor 2 (Supplementary Figure  2A Figure 2B).

Locus 18.
Genotype at rs6516886 correlated with methylation at 6 CpG sites. These are located at either side of the SNP, with cg00065302 (27.4 kb from the SNP), cg05468028 (2.3 kb from the SNP), cg18001427 (1.9 kb from the SNP), and cg20220242 (1.5 kb from the SNP) located on the centromeric side of rs6516886, while cg24751378 (2.7 kb from the SNP) and cg16140273 (62 kb from the SNP) are located on its telomeric side (Supplementary Figure 3A, available on the Arthritis & Rheumatology web site at http://onlin elibr ary.wiley.com/doi/10.1002/ art.40849/ abstract). The CpG site cg00065302 is intergenic and located between LTN1, which codes for E3 ubiquitin-protein ligase listerin, and RWDD2B, which codes for RWD domain-containing protein 2B. The CpG site cg05468028 is located in the gene body of RWDD2, while cg18001427 and cg20220242 are located immediately upstream of this gene. The CpG site cg24751378 is intergenic and located between and upstream of RWDD2B and USP16, which codes for ubiquitin carboxyl-terminal hydrolase 16. The CpG site cg16140273 resides within an intron of MAP-3K7CL, which codes for MAP3K7 C-terminal-like protein. The OA risk-conferring T allele of rs6516886 correlated with increased methylation of cg00065302, cg05468028, cg18001427, cg20220242, and cg24751378, and with decreased methylation of cg16140273 (Supplementary Figure 3B).
Focus on locus 11. Based on the large number of highly significant CpG sites from within locus 11, which all cluster in an interval of only 7.7 kb, we focused our attention on rs11780978 and set out to further investigate this association signal. As noted above, these CpG sites reside within PLEC at chromosome 8q24.3. This region of the genome is gene rich. By analyzing the GTEx database, we first assessed whether there were any known rs11780978 eQTLs operating on genes within the locus. This was the case in a broad range of cell types for PLEC and also for PARP10, GRINA, and SPATC1, which are located in the immediately telomeric region of PLEC ( Figure 1B). The GTEx database does not include data on cartilage tissue samples.    Figures 4B-D).

Methylation QTL, eQTL, and methylationexpression QTL analysis in new patients.
We assessed whether cartilage expression of PLEC, PARP10, or GRINA correlated with rs11780978 genotype and, if so, whether any such eQTLs correlated with methylation. As a first step, we recruited additional OA patients undergoing arthroplasty and extracted DNA and RNA samples from their cartilage. Using pyrosequencing, we measured methylation levels at the PLEC CpG sites in the DNA samples. We focused on 2 of the most significant CpG sites: cg19405177 from the cluster of 3 (cluster 1) and cg14598846 from the cluster of 6 (cluster 2) ( Table 2 and Figure 1). Each pyrosequencing assay captured the relevant CpG site plus additional CpG sites located nearby: 7 sites for cg19405177 and 3 sites for cg14598846. Stratification by rs11780978 genotype replicated our initial mQTL result, and the OA risk-conferring A allele of rs11780978 correlated with higher methylation at cg19405177 (and its 7 flanking CpG sites) and with lower methylation at cg14598846 (and its 3 downstream CpG sites) (Figure 2). We next tested for eQTLs by undertaking AEI analysis in the new patients. For each gene, we selected a transcript SNP that was in linkage disequilibrium with rs11780978, and we identified those rs11780978 heterozygotes that were also heterozygous for 1 or more of the transcript SNPs. AEI was then performed using cDNA synthesized from the cartilage RNA samples. For PLEC we used the transcript SNP rs11783799 (r 2 = 0.93) and analyzed 19 compound heterozygotes. For PARP10 we used the transcript SNP rs11136345 (r 2 = 0.85) and analyzed 22 compound het erozygotes. For GRINA we 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 2 ) 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. used the transcript SNP rs9100 (r 2 = 0.78) and analyzed 11 compound heterozygotes. Significant AEI was observed at PLEC and GRINA, but not at PARP10 (Figure 3). For PLEC, the G allele of rs11783799, which correlates with the OA risk-conferring A allele of rs11780978, demonstrated reduced expression in all 19 compound heterozygotes combined (G:A ratio <1.0; P = 0.02). For GRINA, the T allele of rs9100, which correlates with the A allele of rs11780978, demonstrated increased expression in all 11 compound heterozygotes combined (T:G ratio >1.0; P = 0.005).
We plotted the PLEC and GRINA AEI data for the compound heterozygotes against their individual methylation values. We focused on the PLEC heterozygotes for which we had matched methylation data (up to 15 of the 19 heterozygotes). Additionally, we focused on the GRINA heterozygotes for which we also had methylation data (up to 10 of the 11). For PLEC, 3 cluster-1 CpG sites demonstrated significant correlation between methylation and AEI: CpG1 (P = 0.006), CpG3 (P = 0.031), and CpG6 (P = 0.047) ( Figure 4A). Additionally, we identified a significant correlation between methylation and PLEC AEI at a single clus- 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.  GRINA (B). A, PLEC log 2 allelic expression imbalance (AEI) ratios and B, GRINA log 2 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). ter-2 CpG site (CpG1; P = 0.048). For GRINA, there were no correlations with cluster-1 CpG sites, but 3 cluster-2 CpG sites demonstrated significant correlations: cg14598846 (P = 0.014), CpG2 (P = 0.012), and CpG4 (P = 0.038) ( Figure 4B).
In summary, this analysis replicated the mQTLs, identified cartilage eQTLs at PLEC and GRINA, and revealed the existence of cartilage methylation-expression QTLs (meQTLs) operating on both of these genes but correlating with distinct CpG sites. Methylation at both cluster 1 and 2 was associated with PLEC expression, while only methylation at cluster 2 correlated with GRINA expression.
A genetic difference between clusters 1 and 2. We next sought to assess whether there was a genetic basis for this difference in activity between cluster 1 and cluster 2 CpG sites. Using our original 87 patients, we determined whether any SNPs on the array that were part of the same linkage disequilibrium block as the association signal rs11780978 showed a more significant genotype-methylation correlation with cg19405177 or cg14598846. For cg19405177, we identified rs7819099 (r 2 = 0.83, uncorrected P = 1.66 × 10 −28 versus 2.04 × 10 −20 for rs11780978), and for cg14598846 we identified rs11136336 (r 2 = 0.94, uncorrected P = 3.69 × 10 −27 versus 1.11 × 10 −22 for rs11780978). Both SNPs are also located within PLEC and have a pairwise r 2 value of 0.87. We genotyped these 2 SNPs in our new patients (n = 104) and, as in the 87 patients, we correlated genotype with methylation. In Supplementary Figure 5 (available on the Arthritis & Rheumatology web site at http://onlin elibr ary.wiley.com/doi/10.1002/ art.40849/ abstract), the data are presented as a heatmap for the 87 array patients and for the 104 new patients. For cg19405177, rs7819099 also showed a stronger effect than rs11780978 in the new patients, and in the 2 patient groups combined. For cg14598846, rs11136336 and rs11780978 were more comparable in the new patients and in the 2 groups combined. Overall, these data reveal different polymorphisms as drivers of the methylation-genotype correlations observed at the 2 clusters: rs7819099 for cluster 1 and rs11136336/rs11780978 for cluster 2.
Interaction between the CpG clusters and nearby genes. Finally, we used in silico data and the WashU epigenome browser to search for chromatin interactions between the 2 clusters and nearby genes. No cartilage or chondrocyte data were available for the locus, but we identified relevant data for the human breast cancer cell line MCF7, the chronic myeloid leukemia cell line K562, and the human lymphoblastoid cell line GM12878. For each of these cell lines there was evidence of interaction between the clusters and upstream sequences of several PLEC isoforms (Supplementary Figure 6, available on the Arthritis & Rheumatology web site at http://onlin elibr ary.wiley. com/doi/10.1002/art.40849/ abstract). These data indicate that, at least in these cell lines, these regions harboring the CpG sites physically interact with the PLEC gene.

DISCUSSION
Using genome-wide data, we identified 4 additional OA susceptibility loci in which the risk-conferring allele correlates with differential DNA methylation relative to the nonrisk allele in cis. Combined with our previous analyses (10,12), we have now identified 8 such mQTL loci out of a total of 30 OA association signals studied, which is a frequency of >25%. We suspect that this is an underestimation, considering the low CpG coverage of the Illumina methylation array used (~1.5% of the total CpG sites in the human genome) and our relatively modest sample size of 87 patients. Our analysis, therefore, emphasizes the potentially pivotal mechanistic role that DNA methylation plays in the activity of many OA genetic risk loci, and now warrants targeted editing of the epigenome in order to determine causality. Furthermore, it highlights that this effect is operating on cartilage, the tissue central to the disease process, and is detectable at a clinically relevant point for the patient, namely, when arthroplasty is required.
Of the 4 novel mQTL loci that we report here as significant following FDR correction, we chose to focus on locus 11 because of the large number of highly significant CpG sites at this signal, i.e., 5 of the 9 positive CpG sites had FDR P values of <1 × 10 −10 . A subsequent analysis revealed eQTLs operating on genes within the locus. There is evidence of eQTLs correlating with the OA association SNPs at the other 3 loci (data not shown), and these signals, therefore, also merit follow-up studies.
Methylation QTLs operating at the locus 11 CpG sites in association with the OA SNP, or other variants in high linkage disequilibrium with it, have been previously reported in a range of tissues including adipose, pancreas, lung, lymphocytes, and fibroblasts (34)(35)(36)(37)(38)(39). However, this is the first study of these mQTLs operating in cartilage, and the first study to identify a link between these CpG sites and OA genetic risk. Hypothesizing that the mQTL at this locus signified an eQTL of functional relevance to the association signal, we used a combination of in silico and in vitro analyses to highlight PLEC and GRINA as likely target genes of the rs11780978 association signal.
PLEC codes for plectin, a multifunctional cytoskeletal linker protein that directly interacts with a broad range of cytoplasmic, membrane, and organelle proteins (40,41). Plectin has a range of functions, including acting as a scaffold for signaling proteins and as an organizer of cytoskeletal filaments. The protein is particularly abundant in tissues that are subjected to mechanical stress, and most research to date focuses on plectin in skin, muscle, and blood vessels. Mutations of PLEC result in skin blistering and muscular dystrophy, with the tissue affected and the severity of the disease being determined by the site of the mutation (41). GRINA is a member of a family of genes coding for 6 proteins named "transmembrane BAX inhibitor motif-containing," (TMBIM). These genes encode calcium channels that are present in the Golgi complex, ER, and mitochondria, where they control calcium homeostasis (42). By controlling calcium flow, the TMBIMs | 1295 regulate cell death, including during ER stress, with the majority of the proteins being antiapoptotic (43). GRINA encodes TMBIM3, which is expressed in most cells, and localizes to the Golgi complex and ER, where it suppresses ER stress-induced apoptosis. To the best of our knowledge, plectin and TMBIM3 have not been functionally linked to an arthritis phenotype before.
Our AEI data indicate that the OA risk-conferring A allele of rs11780978 correlates with reduced expression of PLEC but with increased expression of GRINA, while our RNA-Seq data revealed increased PLEC expression but decreased GRINA expression, in OA versus non-OA (femoral neck fracture) cartilage. For PLEC, our interpretation of this is that in an OA joint, additional plectin is required to mitigate the effect of an altered biomechanical load. However, inheritance of the A allele, with its concomitant lower expression of PLEC, attenuates this mitigation. For GRINA, our interpretation is that an alteration in the activity of cells in the OA joint necessitates a reduction in GRINA expression to facilitate controlled cell death. However, the A allele, with its concomitant increased GRINA expression, hinders this response.
One of the most striking outcomes of our investigation was the discovery of 2 clusters of CpG sites that were physically close but had quite separate and distinct characteristics as follows: 1) for cluster 1 the risk-conferring A allele of rs11780978 correlated with higher methylation, whereas for cluster 2 it correlated with lower methylation; 2) CpG sites in both cluster 1 and cluster 2 correlated with an meQTL operating on PLEC, whereas only cluster-2 CpG sites correlated with an meQTL operating on GRINA; and 3) methylation at cluster-1 CpG sites correlated more strongly with genotype at rs7819099, whereas methylation at cluster-2 CpG sites correlated more strongly with genotype at rs11136336/rs11780978. Combined, these data support the notion that at least 2 functional effects are encoded by polymorphisms at the PLEC locus, mediated by methylation on distinct CpG sites and impacting 2 separate genes.
In conclusion, our analysis has identified novel mQTLs and highlighted the interplay between OA genetic risk, DNA methylation, and gene expression in cartilage. We have discovered relevant functional effects on 2 genes from within a single locus and, as such, have elevated these to prime targets for the encoded risk at this locus. These genes and their encoded proteins now merit much more detailed investigation in the context of OA etiology.