Annotating Transcriptional Effects of Genetic Variants in Disease-Relevant Tissue: Transcriptome-Wide Allelic Imbalance in Osteoarthritic Cartilage.

OBJECTIVE: Multiple single-nucleotide polymorphisms (SNPs) conferring susceptibility to osteoarthritis (OA) mark imbalanced expression of positional genes in articular cartilage, reflected by unequally expressed alleles among heterozygotes (allelic imbalance [AI]). We undertook this study to explore the articular cartilage transcriptome from OA patients for AI events to identify putative disease-driving genetic variation.
METHODS: AI was assessed in 42 preserved and 5 lesioned OA cartilage samples (from the Research Arthritis and Articular Cartilage study) for which RNA sequencing data were available. The count fraction of the alternative alleles among the alternative and reference alleles together (φ) was determined for heterozygous individuals. A meta-analysis was performed to generate a meta-φ and P value for each SNP with a false discovery rate (FDR) correction for multiple comparisons. To further validate AI events, we explored them as a function of multiple additional OA features.
RESULTS: We observed a total of 2,070 SNPs that consistently marked AI of 1,031 unique genes in articular cartilage. Of these genes, 49 were found to be significantly differentially expressed (fold change <0.5 or >2, FDR <0.05) between preserved and paired lesioned cartilage, and 18 had previously been reported to confer susceptibility to OA and/or related phenotypes. Moreover, we identified notable highly significant AI SNPs in the CRLF1, WWP2, and RPS3 genes that were related to multiple OA features.
CONCLUSION: We present a framework and resulting data set for researchers in the OA research field to probe for disease-relevant genetic variation that affects gene expression in pivotal disease-affected tissue. This likely includes putative novel compelling OA risk genes such as CRLF1, WWP2, and RPS3.

Introduction

Due to the increased proportion of elderly persons in the human population, osteoarthritis (OA) has become one of the major musculoskeletal diseases 1. While all joint tissues have been implicated in OA pathology, the disease is characterized primarily by progressive degradation and calcification of articular cartilage 2. Both gene‐targeted research 3, 4, 5 and genome‐wide research 6, 7, 8, 9 showed that a multitude of genes are involved in the currently irreversible destruction of articular cartilage that precedes total joint replacement surgery, which is at present the only effective treatment for end‐stage OA. In this regard, numerous studies have shown altered regulation of gene expression that reflects, attenuates, and/or stimulates OA‐mediated cartilage degradation 10, 11, 12, 13. Moreover, multiple OA risk alleles of single‐nucleotide polymorphisms (SNPs) were shown to consistently modulate OA pathology by altering transcription of the respective genes in articular cartilage, commonly referred to as an allelic imbalance (AI) 14, 15, 16, 17, 18, 19. Notable recent examples are the genes ALDH1A2 18 and MGP 20. Hence, it is clear that in‐cis genetic regulation of transcription plays a substantial role in cartilage homeostasis and, therefore, in OA pathophysiology.

Despite the evidence for in‐cis genetic regulation of transcription in OA susceptibility, genome‐wide association studies (GWAS) have thus far failed to explain the larger part of the hereditary component of OA 21. In this regard, a large number of the tested SNPs in GWAS likely bear no biologic function in relation to the addressed phenotype or disease‐relevant tissues 22, resulting in massive inflation of possibly biologically irrelevant statistical tests and thus the multiple testing correction penalty. Consequently, large numbers of SNPs that do bear biologic functionality in the context of OA are missed. Furthermore, SNPs that reside within linkage disequilibrium (LD) blocks are hard to interpret, as association analysis is inherently unable to distinguish disease‐relevant alleles from merely statistically associated alleles.

In previous studies, we and others have used targeted approaches to address AI events of putative as well as established OA susceptibility genes 14, 16, 17, 18, 23, 24. Given the successful identification of the transcriptional consequences of multiple OA‐associated SNPs, we have aimed to characterize, on a transcriptome‐wide scale, novel SNPs that tag AI of genes expressed in articular cartilage, and we have subsequently identified those that appear to confer susceptibility to OA. Finally, further using the RNA sequencing data set, we ran analyses to identify AI genes whose expression was additionally modified with severity of OA pathophysiology as reflected by differential expression between preserved and lesioned OA cartilage.

Materials and Methods

Cohorts

Ethical approval for the Research Arthritis and Articular Cartilage (RAAK) study was obtained from the medical ethics committee of the Leiden University Medical Center (P08.239), and informed consent was obtained from all patients included. For the current study, RNA sequencing data were available from preserved and lesioned cartilage from 21 patients (6 with hip OA and 15 with knee OA), complemented by an additional 21 preserved samples (from 14 patients with hip OA and 7 patients with knee OA) and 5 lesioned samples (from 2 patients with hip OA and 3 patients with knee OA) (see Supplementary Table 1, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract). For cartilage sampling details, see refs. 19 and 25.

RNA sequencing data

After RNA isolation (RNeasy Mini Kit, RNA integrity number >7; Qiagen), paired‐end 100‐bp RNA library sequencing (Illumina TruSeq RNA Library Prep Kit, Illumina HiSeq 2000) resulted in an average of 10 million clusters. Reads were aligned using GSNAP (an R package within Bioconductor; https://rdrr.io/bioc/gmapR/) against the human (hg19) reference genome, while known Dutch SNPs (Genome of the Netherlands) were masked to aid in preventing potential reference alignment bias. AI events were assessed on SNPs called using SNVMix2 with default settings 26 with minimum coverage of 25 and at least 10 reads (R) per allele. AI is reported as the average fraction (φ) of the alternative allele reads (R ) among the total number of reads (R = R + R ) at the position of the respective genetic variation per sample (i):To detect SNPs that robustly mark imbalance, 2 binomial tests were performed per heterozygote and per SNP under the null hypothesis that the amount of imbalance is either greater or smaller than 0.477. A meta‐analysis (meta; http://www.r-project.org/) per SNP among heterozygous individuals (null hypothesis median φ = 0.49) was performed to generate a meta‐φ and P value per SNP with a false discovery rate (FDR) correction for multiple comparisons. To allow independent samples in the meta‐analysis only, the preserved cartilage of each sample pair (n = 21) was used complemented with the individual preserved (n = 21) and lesioned (n = 5) OA samples.

Unfortunately, SNVMix2 discards strand specificity. While in general this does not pose an issue when annotating the AI direction to an effector allele, it does so for A>T, T>A, G>C, and C>G SNPs. Therefore, for these SNPs (n = 119) we supplied (see Supplementary Table 2, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract) minor allele frequencies (MAFs) and SNPs in strong LD for GWAS look‐ups. Using the edgeR package, fragments per gene were used to assess the dispersion by quantile‐adjusted conditional maximum likelihood 27. Subsequently, differential gene expression analysis was performed pairwise between preserved and lesioned samples for which we had RNA of both (n = 21) (Table 1) followed by FDR correction. Gene Ontology (GO) term enrichment analysis was performed using the tool DAVID available online 28.

Table 1

Sample characteristics of preserved and lesioned OA articular cartilage in the Research Arthritis and Articular Cartilage studya

Tissue type	No. of samples	Age, mean ± SD years	No. of men	No. of women	OA articular cartilage
					Preserved	Lesioned	Preserved–lesioned pairs
Knee	25	69 ± 9	4	21	22	18	15
Hip	22	66 ± 9	5	17	20	8	6
All	47	68 ± 9	9	38	42	26	21

OA = osteoarthritic.

Genotype data

Using Illumina HumanOmniExpressExome chips, genome‐wide genotyping data were constructed for 216 samples from the RAAK study. SNPs with <95% call rate, Hardy‐Weinberg equilibrium <10−4, MAF <0.01, or located on the sex chromosomes were removed prior to imputation against the 1000 GenomesV3 March 2012 reference panel 29. We removed SNPs for which the imputation quality of 0.4 was not met 12, 30.

TaqMan assay

Conventional TaqMan genotyping was performed on both genomic DNA and articular cartilage complementary DNA (cDNA) 31 from 6 patients (2 females and 4 males) who underwent total joint replacement surgery of the knee due to primary OA. An allele‐specific custom TaqMan assay for rs7256319 (ThermoFisher Scientific) was used to quantify the allele ratio in cDNA samples. AI of cDNA was normalized against the genomic DNA ratio (with an inherent 1:1 allele ratio) as a reference.

Results

Transcriptome‐wide discovery of articular cartilage AI events

To understand how genetic variation contributes in‐cis to transcriptional regulation in articular cartilage on a transcriptome‐wide scale, we first called heterozygous SNPs (dbSNP144) using RNA sequencing data from articular cartilage derived from patients who underwent total replacement surgery of either the hip joint (n = 22) or knee joint (n = 25) due to primary OA (Table 1; also see Supplementary Table 1, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract). After filtering by the number of read counts per position (R ≥10, R ≥10, and R ≥25), selecting for heterozygous SNPs present in at least 2 individuals, removing SNPs present in multiple or no distinct transcripts, and discarding the HLA locus, we defined φ for 13,853 SNPs as the measure of imbalance (Figure 1A), which denotes the fraction of R among R . Possibly due to reference bias, a considerable number of SNPs were marked as allelic imbalanced by φ <0.1 or φ >0.9 (n = 418) and were subsequently removed prior to further analyses. As such, in Supplementary Table 2 (http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract) we show AI, defined as φ across all individuals heterozygous for the 13,435 SNPs.

Figure 1

Distribution of allelic imbalance (AI) events in articular cartilage. AI is reported as the average fraction (φ) of the alternative allele reads among the total number of reads. A, AI was defined for 13,853 called variants after selecting for at least 2 heterozygotes, selecting single‐nucleotide polymorphisms (SNPs) present in only single genes and removing low counts. B, After filtering by allelic fraction (0.1<φ<0.9) and meta false discovery rate <0.05, 2,070 SNPs remained that marked significant AI of 1,031 unique genes.

Subsequently, a meta‐analysis among heterozygous individuals of each SNP (null hypothesis: median φ = 0.49) and subsequent FDR correction for multiple testing revealed 2,070 SNPs that significantly marked AI among 1,031 genes (Figure 1B; also see Supplementary Table 3, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract). To allow unambiguous annotation of the AI direction to a GWAS effector allele, for A>T, T>A, G>C, and C>G SNPs with MAFs surrounding 0.50 (n = 119), we supplied MAFs and SNPs in strong LD.

Intersection of genes subject to AI with those differentially expressed between preserved and paired OA‐lesioned cartilage

While articular cartilage genes subject to AI due to genetic variation could contribute to OA pathophysiology in various ways (e.g., in cartilage development or homeostasis), it can be expected that allelic imbalanced genes that additionally mark the articular cartilage’s disease state are more likely to contribute to or attenuate disease progression. Therefore, we went back to the original expression data and determined differential expression in patients for whom paired RNA sequencing data of both preserved and OA‐lesioned articular cartilage were generated (6 hip joints and 15 knee joints) (Table 1).

Of the 10,468 Ensembl gene identifiers with at least 5 counts per million, 137 and 86 were observed to be significantly (FDR <0.05) down‐regulated (fold change <0.5) and up‐regulated (fold change >2), respectively, in lesioned cartilage compared to preserved cartilage (see Supplementary Figure 1A and Supplementary Table 4, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract). As has been shown by microarray studies that have used a similar design, subsequent GO term enrichment analysis (see Supplementary Table 5, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract) revealed significant enrichment for inflammatory pathways (e.g., SCUBE1, CFH, and CXCL14) (see Supplementary Figures 1B–D, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract), pathways of response to wound healing (e.g., NOTCH3, BMP5, and SERPINE1) (see Supplementary Figures 1E–G), and joint development–associated pathways (e.g., SPP1, MMP3, and COL9A1) (see Supplementary Figures 1H–J).

Of the 223 differentially expressed genes, 49 were additionally subject to AI, marked by 128 SNPs (see Supplementary Table 6, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract). A notable example was the CRLF1 gene, which was subject to highly consistent AI, with the T allele of rs7256319 marking consistently lower expression of CRLF1 compared to the reference allele C (φ = 0.29, FDR = 4.02 × 10−21) (Figure 2A). Moreover, as shown in Figure 2B, the AI of CRLF1 was confirmed by custom TaqMan assay performed in 5 preserved and 5 lesioned articular cartilage samples, originating from 6 independent patients who underwent total knee replacement surgery. In parallel, expression of CRLF1 in the current data set differed significantly between preserved and OA‐lesioned cartilage, with significant up‐regulation in OA‐affected cartilage (fold change 4.6, FDR = 3.1 × 10−10) (Figure 2C).

Figure 2

expression is subject to allelic imbalance (AI) and is modulated in articular cartilage with osteoarthritis (OA)–induced destruction. A, AI of rs7256319 in RNA sequencing data set, reported as the average fraction (φ) of the alternative allele reads among the total number of reads with 95% confidence interval, showing decreasing expression of transcript of the alternative T allele in preserved and lesioned OA cartilage relative to the reference C allele (φ = 0.29, false discovery rate [FDR] = 4.02 × 10−21). Preserved OA cartilage is depicted in blue, lesioned OA cartilage in red, and meta‐φ (META) in black. B, Replication of AI expression with rs7256319 by TaqMan genotyping in 6 additional knee samples, confirming the observed lower expression of the alternative allele T relative to the reference allele C of rs7256319 (T:C ratio = 0.63). Preserved OA cartilage is depicted in blue, lesioned OA cartilage in red, and genomic DNA (used as the reference ratio) in black. The meta‐φ of genomic DNA, lesioned cartilage complementary DNA (cDNA), and preserved cartilage cDNA is depicted in black, red, and blue, respectively. Data are shown as box plots. Each box represents the 25th to 75th percentiles. Lines inside the boxes represent the median. Lines outside the boxes represent the 10th and 90th percentiles. In A and B, horizontal dashed lines depict equal ratios of the CRLF1 rs7256319 alleles. C, Differential expression analyses of showing significantly up‐regulated expression in lesioned (red) compared to paired preserved (blue) OA articular cartilage (fold change 4.6, FDR‐corrected P = 3.1 × 10−10). RAAK = Research Arthritis and Articular Cartilage (study).

Cartilage AI SNPs that contribute to OA susceptibility

On the basis of the significant AI SNPs in articular cartilage (see Supplementary Table 3, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract), it is to be expected that these SNPs are enriched for those conferring genetic susceptibility to OA. Hence, for 173 AI SNPs with FDR <5 × 10−8, not previously reported as OA risk genes, we went to the genome‐wide association catalogs of OA to look for their genetic association signal. We used the GWAS meta‐analysis for hip OA performed under the auspices of the Translational Research in Europe Applied Technologies in Osteoarthritis (TreatOA) consortium 32 and the GWAS meta‐analysis on cartilage thickness (measured by the OA endophenotype minimal joint space width [JSW]) 8. Given that the entire genome was not being assessed for genetic association, a nominal genetic association (P < 0.05) was considered. Moreover, we only checked for genetic association when the AI SNP could be used directly or when a clear proxy SNP was identified.

As shown in Table 2, we observed multiple AI SNPs conferring susceptibility to hip OA and/or minimal JSW. Notable examples in Table 2 are the AI SNP rs3133187 in the RPS3 gene conferring the most significant susceptibility to hip OA (Figure 3A) and the AI SNP rs1052429 in the WWP2 gene associated with minimal JSW (Figure 3B). Alternatively, in Supplementary Tables [Link], [Link] (http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract), we show the intersection of the 65 AI SNPs significantly associated with hip OA in the TreatOA GWAS meta‐analysis (FDR <0.05) and the 63 AI SNPs significantly associated with minimal JSW in the GWAS meta‐analysis on cartilage thickness (P < 0.05). Among these lists we find SNPs in compelling known (MGP, ALDH1A2, FRZB, COL11A1, PLEC) and potentially novel (ACAN, MATN, TNC, VEGFA, PLOD2) OA risk genes.

Table 2

Intersection of highly significant allelic imbalance SNPs with published genome‐wide association data sets of minimal JSW and hip OAa

Allelic imbalance						Genome‐wide association
SNP	Positional gene			Meta‐φ	FDR	SNP	R2	Minimal JSW					Hip OA
		AA	AAF					EA	EAF	Beta	SEM	P	EA	EAF	Beta	SEM	P
rs4744	PLA2G2A	A	0.11	0.46	8.1 × 10−134	rs11677	1.0	A	0.11	−0.05	0.02	0.0091	A	0.12	0.00	0.04	0.9360
rs4605	FMOD	C	0.70	0.49	0.000	rs4605	NA	C	0.52	−0.03	0.01	0.0679	C	0.46	−0.09	0.03	0.0057
rs3190	PPP1CB	G	0.56	0.59	6.3 × 10−10	rs3190	NA	G	0.57	0.03	0.01	0.0076	G	0.57	−0.00	0.02	0.8820
rs1054629	IBSP	T	0.18	0.41	1.0 × 10−24	rs13144371	1.0	A	0.24	0.03	0.01	0.0340	A	0.24	−0.05	0.03	0.0801
rs6546	C1QTNF3	G	0.30	0.39	4.8 × 10−13	rs840385	0.88	C	0.40	−0.03	0.01	0.0168	C	0.39	−0.03	0.03	0.2361
rs3549	SPARC	C	0.39	0.45	1.1 × 10−58	rs1059829	1.0	A	0.53	0.03	0.01	0.0284	A	0.53	−0.02	0.02	0.3896
rs3829078	CA9	G	0.05	0.59	2.5 × 10−14	rs3829078	NA	G	0.11	0.04	0.02	0.0495	G	0.11	−0.04	0.04	0.2710
rs13321	TNC	G	0.67	0.53	2.4 × 10−16	rs12347433	0.1	T	0.72	−0.02	0.01	0.1487	T	0.72	0.06	0.03	0.0256
rs1871452	CHST3	A	0.60	0.42	5.4 × 10−25	rs731027	1.0	T	0.52	0.01	0.01	0.3073	T	0.50	−0.06	0.02	0.0087
rs3133187	RPS3	G	0.12	0.42	5.0 × 10−9	rs3133187	NA	G	0.06	−0.04	0.03	0.1420	G	0.06	0.18	0.05	0.0011
rs1800801	MGP	T	0.74	0.41	2.6 × 10−26	rs4236	0.9	T	0.61	−0.03	0.01	0.0433	T	0.60	0.03	0.03	0.2068
rs3737548	COL2A1	T	0.23	0.53	5.1 × 10−72	rs1635553	0.4	G	0.46	−0.02	0.01	0.2290	G	0.46	−0.06	0.03	0.0161
rs6647	SERPINA1	G	0.75	0.41	7.2 × 10−283	rs6647	NA	G	0.77	−0.03	0.01	0.0848	A	0.76	0.06	0.03	0.0388
rs1052429	WWP2	A	0.86	0.56	3.4 × 10−55	rs1566452	0.9	A	0.73	−0.04	0.01	0.0028	A	0.71	−0.05	0.03	0.0656
rs2646108	CRISPLD2	A	0.18	0.63	2.6 × 10−9	rs2646108	NA	A	0.20	−0.01	0.02	0.7410	A	0.19	0.09	0.03	0.0053
rs6554	UBA52	T	0.70	0.43	1.1 × 10−53	rs6554	NA	T	0.60	0.00	0.01	0.9300	T	0.61	0.05	0.03	0.0416
rs7499	COL18A1	A	0.47	0.63	2.0 × 10−11	rs7499	NA	A	0.40	−0.01	0.01	0.3528	A	0.40	−0.05	0.03	0.0332
rs3088026	COL6A2	T	0.07	0.53	5.5 × 10−3	rs3088026	NA	T	0.09	−0.02	0.02	0.4125	T	0.09	−0.12	0.04	0.0055

SNPs = single‐nucleotide polymorphisms; JSW = joint space width; OA = osteoarthritis; AA = alternative allele; AAF = alternative allele frequency; φ = count fraction of the alternative allele reads among the total number of reads; FDR = false discovery rate relative to median φ = 0.49 as reference; R2 = measure of linkage disequilibrium of proxy SNP with allelic imbalance SNP; EA = effect allele; EAF = effect allele frequency; NA = not applicable.

Figure 3

Compelling cartilage‐specific allelic imbalance (AI) genes. AI is reported as the average fraction (φ) of the alternative allele reads among the total number of reads with 95% confidence interval. A, Significant (P = 4.97 × 10−9) allelic imbalanced expression at rs3133187 with the G allele decreasing expression of transcript in preserved and lesioned osteoarthritic (OA) cartilage relative to the A allele (meta‐φ = 0.42). The G allele of rs3133187 additionally confers significant association with hip OA (P = 1.1 × 10−3). B, Significant (P = 3.37 × 10−55) allelic imbalanced expression at rs1052429 with the A allele increasing expression of transcript in preserved and lesioned OA cartilage relative to the G allele (meta‐φ = 0.56). The A allele of rs1052429, linked to the A allele of rs1566452 (r2 = 0.92, D′ = 1.0), additionally confers significant association with reduced minimal joint space width (P = 2.8 × 10−3), resulting in cartilage degeneration or OA. Horizontal dashed lines depict equal ratios of alleles. RAAK = Research Arthritis and Articular Cartilage (study).

Cartilage AI SNPs that contribute to expression quantitative trait loci (eQTLs)

To provide additional supporting data on the AI findings in the putative OA risk SNPs in Table 2, we combined the RNA sequencing data of the genes with genome‐wide SNP data (Illumina HumanOmniExpressExome) of RAAK study samples to extract gene‐targeted cartilage eQTLs on the basis of 50 samples that were overlapping. Additionally, we explored eQTL data of the Genotype‐Tissue Expression (GTEx) Project. The fact that GTEx data are merely from tissues other than those particularly relevant to OA indicates the generalizability of the identified AI to other tissues. As shown in Supplementary Table 9 (http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract), the GTEx eQTL data largely support (the direction of) the identified AI effects also in other tissues.

OA susceptibility SNPs that show AI in cartilage

Since identified OA risk SNPs have been demonstrated to frequently confer risk by modifying expression of positional genes in‐cis, the aforementioned genome‐wide AI data set of articular cartilage (see Supplementary Table 3, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract) can function as a database to assess in silico the direction of effect of such identified susceptibility SNPs. Hence, we reviewed the literature to intersect the list of currently published robust genetic OA association signals with the identified AI SNPs. We included robustly identified susceptibility SNPs and their positional genes associated with knee, hand, and hip OA as well as OA‐associated phenotypes such as minimal JSW and markers of cartilage turnover (soluble cartilage oligomeric matrix protein). Of the 46 genes linked to the reported OA risk SNPs, 36 genes were detected in our data set (i.e., had detectable expression levels in cartilage) and were heterozygous carriers of a coding SNP. As shown in Supplementary Table 10 (http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract), we found 11 previously identified OA risk genes to contain SNPs that mark AI in articular cartilage. Of these, the respective risk alleles marked lower expression through AI in heterozygotes of FRZB (rs7775), COL11A1 (rs2615977 and rs1676486), IGFBP3 (rs788748 through double heterozygotes with rs6670), ALDH1A2 (rs3204689), and MGP (rs4764133 in LD with rs1800801). For SMAD3, CDCL5, COL12A1, BCAP29, PIK3R1, and COMP, AI was detected but not in relation to the reported risk alleles.

To provide additional functional data to these previously identified OA risk genes, level and differential expression data between preserved and lesioned OA articular cartilage are depicted in Supplementary Table 11 (http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract). For IGFBP3, PIK3R1, BCAP29, COL12A1, FRZB, ALDH1A2, and MGP, significant differential expression (FDR <0.05) was observed in addition to AI. Finally, to find supporting data on the AI findings in the OA risk SNPs, we combined the RNA sequencing data of the OA risk genes with genome‐wide association data (Illumina HumanOmniExpressExome) of RAAK study samples to extract gene‐targeted eQTL data on the basis of 50 cartilage samples that were overlapping. Additionally, to the eQTL data in Supplementary Table 11, we added data from the GTEx Project, although from tissues other than those particularly relevant to OA. For the OA genes GNL3, FTO, NCOA3, MICAL3, IFRD1, IGFBP3, TGFA, and MGP, we found significant eQTL effects (P < 0.05) of risk SNPs that substantiated the respective AI effects.

Discussion

Our approach in the current study comprises a concept framework for complex traits to identify disease‐relevant genetic variation, as reflected by allele‐associated transcription levels in cartilage, a pivotal tissue in the disease. We have aimed to present the reported observations as a legacy data set for researchers in the field to probe for their gene or SNP of interest. Herein, we highlight notable examples. Among our (highly) significant AI SNPs we confirmed well‐known genes that have previously been reported by others to confer robust risk of OA (e.g., MGP, ALDH1A2, IGFBP3, and FRZB) 17, 33, 34. We hypothesize that among our (highly) significant AI SNPs and particularly those that show additional, differential expression between preserved and lesioned OA cartilage (e.g., CRLF1 [Figure 2; also see Supplementary Table 6, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract]) or genetic association with OA phenotypes (e.g., WWP2 and RPS3 [Figure 3 and Table 2; also see Supplementary Tables [Link], [Link], http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract]) are putative novel compelling OA risk genes.

The CRLF1 gene, encoding for cytokine receptor–like factor 1 protein, harbors the rs7256319 C>T SNP that has marked imbalanced expression of its respective alleles in articular cartilage, reflected by consistent lower expression of the alternative allele T in comparison with the reference allele C among heterozygotes. As was also reported previously 25, 35, CRLF1 appeared to be significantly up‐regulated in lesioned compared to preserved OA articular cartilage (fold change 4.6, FDR = 3.1 × 10−10), as was its signaling partner CLCF1 (fold change 2.1, FDR = 1.0 × 10−6), while the protein complex signaling receptor gene CNTFR was significantly down‐regulated (fold change 0.3, FDR = 1.9 × 10−8) (see Supplementary Table 3, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract). Additionally, it was shown by Tsuritani et al 35 that up‐regulation of the cytokine receptor–like factor 1/cardiotrophin‐like cytokine complex in ATDC5 cells disrupts cartilage homeostasis and promotes progression of OA by enhancing the proliferation of chondrocytes and suppressing the production of cartilage matrix. As such, we hypothesize that the alternative allele T of rs7256319 in heterozygote carriers may be able to mitigate CRLF1/CLCF1 signaling toward ongoing cartilage degradation due to primary OA.

Among the notable novel putative OA genes in Table 2 is WWP2, which showed (in addition to multiple coding SNPs marking consistent AI expression of WWP2 in cartilage) significant differential expression between preserved and lesioned OA cartilage (fold change 0.78, FDR = 0.0053; results not shown) and a signal of genetic association with minimal JSW (P = 0.0028) (Figure 3). Based on these data, we hypothesize that allele rs1052429 A located in WWP2 is an OA susceptibility allele that acts via higher expression of WWP2 in cartilage, which is associated with lower minimal JSW and thus with degeneration of cartilage. In addition, expression of WWP2 was previously shown to be consistently modified by methylation at CpG sites 36. The WWP2 protein is a member of the Nedd4 family of E3 ligases, which play an important role in protein ubiquitination 37. Moreover, the encoded protein was shown to physically interact with SOX9, thereby affecting the transcriptional activity of SOX9 via translocation to the nucleus 38.

Similarly, RPS3, encoding ribosomal protein S3, is subject to significant AI marked by rs3133187 with the G allele decreasing expression of RPS3 transcript in cartilage (FDR = 4.97 × 10−9) and significant association with hip OA (P = 1.1 × 10−3). Ribosomal protein S3 is a multifunctional DNA repair endonuclease and ribosomal protein, yet its role in inducing apoptosis through activation of CASP8/CASP3 may be relevant to mention with respect to OA 39. It was recently shown by Jeon et al 40 that selective removal of senescent cells may attenuate the development of OA. We hypothesize that decreasing RPS3 expression by the OA risk allele G of rs3133187 negatively affects removal of senescent chondrocytes by cytokine‐induced apoptosis. Nonetheless, functional studies are necessary to verify the exact mechanism by which the AI in CRLF1, WWP2, and RPS3 contributes to cartilage degeneration in humans. It would be preferable to investigate this using human in vitro micromass cultures in which the expression of the genes (e.g., by lentiviral induction) is modified in the direction of the risk alleles as previously found with, for example, GDF5 15 or DIO2 14.

Our RNA sequencing data set consisted of both preserved and lesioned OA samples, as well as knee and hip cartilage. Due to our focus on significant AI effects across our samples, we have ignored possible variation in the AI effects between preserved and lesioned OA cartilage or between knee and hip joints. Such variation could arise, e.g., due to local expression of transcription factors that exaggerate or attenuate the observed AI. By providing AI, defined by φ across all individuals heterozygous for any of the 13,435 SNPs (see Supplementary Table 2, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract), possible differences in AI for preserved/OA or hip/knee cartilage may be explored. It should be mentioned, however, that the current sample size and the variable number of heterozygous individuals likely precludes a robust statement on either analyses and will require additional targeted AI measurements by, for example, TaqMan assay.

By combining reported genetic OA signals with AI SNPs that alter transcription of articular cartilage genes in ‐ cis, we also had the opportunity to address functionality of OA risk alleles. For example, the A allele of rs788748, located upstream of IGFBP3, is associated with lower odds of hip OA 41. Given that this SNP is not located in an exon, we assessed its potential revelance to AI through rs6670 double heterozygotes, which revealed that the protective rs788748 A allele marks lower expression of IGFBP3 compared to the G allele. In addition, significant up‐regulation of the gene was observed in OA‐affected cartilage compared to paired preserved cartilage, as well as a similar eQTL effect of the rs788748 SNP (see Supplementary Table 9, http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract), similar to what has been reported by others 8. These observations solidify the notion of IGFBP3's role in OA, and we propose that the reported protective effect of the rs788748 A allele is mediated through lower expression of the transcript.

For a selected number of comparable (although gene‐targeted) studies, we discuss the respective confirmations and/or discrepancies. A number of SNPs are known to mark AI in articular cartilage, as has been shown by gene‐targeted approaches. We were able to replicate the earlier observed AI of ALDH1A2 18 and MGP 20 and to a lesser extent that of DIO2 14 and GDF5 42. Furthermore, for rs11177 (GNL3) 16 and rs6617 (SPCS1) 16 in our data set, AI was not as obvious as previously reported, and we did not observe heterozygotes for rs143383 (GDF5) 23 or rs3815148 (HBP1) 24.

For COL11A1, AI was previously thoroughly investigated by Raine et al 17 in view of the OA risk SNP rs2615977 7 and the lumbar disc herniation SNP rs1676486 43. They showed that considerable AI of COL11A1 was correlated with the lumbar disc herniation risk allele of SNP rs1676486, while the observed AI of COL11A1 was not correlated with the OA risk allele of SNP rs2615977. The latter result was based on the AI effect observed in individuals double heterozygous for rs2615977 and rs9659030. As shown in Supplementary Figure 2 and Supplementary Table 8 (http://onlinelibrary.wiley.com/doi/10.1002/art.40748/abstract), we observed 2 independent COL11A1 SNPs (D′ = 0.2, r2 = 0.01) 44 with highly significant meta‐φ AI. The first was the abovementioned SNP rs9659030 with meta‐φ AI of 0.65 (P = 1.2 × 10−25) and the second was SNP rs2229783 with meta‐φ AI of 0.53 (P = 3.7 × 10−8). Notably, and in contrast to the results obtained by Raine et al 17, the extent and consistency of AI for rs9659030 were considerably higher than those for rs2229783 and rs1676486. The suggestive evidence of AI of the lumbar disc herniation risk allele of SNP rs1676486 as well as the AI of its proxy SNP rs2229783 (D′ = 1, r2 = 0.3) confirm that the lumbar disc herniation risk allele is likely associated with lower COL11A1 expression. Moreover, the AI of the noncoding OA risk SNP rs2615977, which is not in strong LD with a coding SNP, remains unclear. Yet, based on the proxy SNP rs1031820 (D′ = 0.77, r2 = 0.16), it may be that the OA risk allele of rs2615977 also acts via lower expression of COL11A1.

Taken together, these confirmations and/or discrepancies indicate first and foremost that additional replication to verify AI is required to increase confidence in the observed AI (e.g., by better preselection on heterozygous samples). Second, it stresses the fact that observed AI reflects regulatory properties of the respective LD block and does not per se identify genetic variation that affects respective gene expression levels mechanistically. Furthermore, despite the applied filtering steps and statistics, the list of significant AI SNPs potentially contains a number of false positives, some of which could have originated from alignment bias. While future novel alignment and other bioinformatic approaches 45 might address these issues from a more fundamental perspective, in the present study we have aimed to reduce false‐positive AI SNPs by including multiple filtering steps (0.1<φ<0.9, include SNPs with at least 2 heterozygotes, and null hypothesis adjustment).

While in canonical GWAS a strict genome‐wide significance level of 5 × 10−8 is imposed due to the vast amount of SNPs that are tested for, we postulate that providing SNPs that are more likely to affect expression of genes in‐cis in a disease‐relevant tissue could aid the search for the functional susceptibility SNPs and the putative OA risk gene. Nevertheless, among the significant AI SNPs, we did not necessarily obtain a clear enrichment of putative significant OA risk SNPs; among the SNPs for which we had both AI and GWAS data, 6.7% showed significant association with the OA phenotypes, while among the significant AI genes (FDR <0.05), 7.2% showed significant association with the OA phenotypes. Further downstream selection criteria, such as (but not limited to) significant differential expression between preserved and OA‐lesioned cartilage and/or trans‐eQTL analysis, will help tailor genetic association analyses even more and might attribute SNPs to specific disease facets, such as extent of cartilage degradation, as we have shown in the present study. Of note, inherent to our study design, we might have missed genes that affect joint morphology or cartilage integrity during development and/or that change expression in healthy cartilage or during early‐stage OA.

In summary, we present a framework and resulting data set for researchers in the OA research field to probe for disease‐relevant genetic variation that affects gene expression in pivotal disease‐affected tissue. This likely includes putative novel compelling OA risk genes such as CRLF1, WWP2, and RPS3.

Author Contributions

All authors were involved in drafting the article or revising it critically for important intellectual content, and all authors approved the final version to be published. Dr. Meulenbelt had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Study conception and design

Den Hollander, Couthino de Almeida, Lakenberg, ‘t Hoen, Ramos, Meulenbelt.

Acquisition of data

Den Hollander, Bomer, Nelissen, Ramos, Meulenbelt.

Analysis and interpretation of data

Den Hollander, Pulyakhina, Boer, van der Breggen, Arindrarto, Couthino de Almeida, Sentner, Laros, ‘t Hoen, Slagboom, van Meurs, Ramos, Meulenbelt.

Supplementary Figures

Click here for additional data file.

Supplementary Table 1

Click here for additional data file.

Supplementary Table 2

Click here for additional data file.

Supplementary Table 3

Click here for additional data file.

Supplementary Table 4

Click here for additional data file.

Supplementary Table 5

Click here for additional data file.

Supplementary Table 6

Click here for additional data file.

Supplementary Table 7

Click here for additional data file.

Supplementary Table 8

Click here for additional data file.

Supplementary Table 9

Click here for additional data file.

Supplementary Table 10

Click here for additional data file.

Supplementary Table 11

Click here for additional data file.