Expression quantitative trait loci analysis in rheumatoid arthritis identifies tissue specific variants associated with severity and outcome

Abstract

Objective: Genome-wide association studies have successfully identified more than 100 loci associated with susceptibility to rheumatoid arthritis (RA). However, our understanding of the functional effects of genetic variants in causing RA and their effects on disease severity and response to treatment remains limited.

Methods: In this study, we conducted expression quantitative trait locus (eQTL) analysis to dissect the link between genetic variants and gene expression comparing the disease tissue against blood using RNA-Sequencing of synovial biopsies (n=85) and blood samples (n=51) from treatment-naïve patients with RA from the Pathobiology of Early Arthritis Cohort.

Results: This identified 898 eQTL genes in synovium and genes loci in blood, with 232 genes in common to both synovium and blood, although notably many eQTL were tissue specific. Examining the HLA region, we uncovered a specific eQTL at HLA-DPB2 with the critical triad of single-nucleotide polymorphisms (SNPs) rs3128921 driving synovial HLA-DPB2 expression, and both rs3128921 and HLA-DPB2 gene expression correlating with clinical severity and increasing probability of the lympho-myeloid pathotype.

Conclusions: This analysis highlights the need to explore functional consequences of genetic associations in disease tissue. HLA-DPB2 SNP rs3128921 could potentially be used to stratify patients to more aggressive treatment immediately at diagnosis.

Keywords: Arthritis, Rheumatoid; Methotrexate; Pharmacogenetics; Polymorphism, Genetic; Synovitis.

Figure 1

cis-eQTL for Blood and Synovium in Early Rheumatoid Arthritis cis-eQTL analysis was performed on 85 synovial and 51 blood samples using matrixEQTL (A). P values were calculated based on the t-statistic from a linear model with four PCA eigenvalues and four PEER factors as covariates. Following filtering of SNPs based on minor allele frequency (maf≥0.05), and genome wide significance (FDR≤0.01), 898 unique genes were found to have associated eQTL in synovium, and 1251 in blood. Log quantile-quantile (QQ) p value plots (B) indicate the overall level of significance from cis-eQTLs in each tissue. Manhattan plots (C) show genome-wide cis-interactions for variants located within ±0.5 Mb of genes in synovium (top) and blood (bottom). The dashed line represents genome-wide significance at FDR=0.01. Genes were labelled where there were over 50 significant gene-SNP interactions. eQTL, expression quantitative trait locus; FDR, false discovery rate; PCA, principal component analysis; PEER, Probabilistic Estimation of Expression Residuals; SNP, single-nucleotide polymorphism.

Figure 2

Selection of eQTL SNPs and genes associated with clinical response. Columns are ordered from left to right. The first column (far left) shows locus Manhattan plots for representative eQTLs with -log₁₀ p value on the y axis and SNP chromosomal position on the x axis. P values were calculated based on the t-statistic from a linear model with four genotype eigenvalues and four RNA-Seq PEER factors as covariates. Synovial eQTL at the PEX6 and ERAP2 loci, and blood eQTL at the C3AR1 and SPG20 loci are shown. The second column shows violin plots of expression of the most strongly associated eQTL SNP and baseline RNA-Seq gene expression in either synovium or blood (statistical analysis by linear model with four PCA andfour PEER factor covariates). The third column shows relationships between baseline RNA-Seq gene expression and clinical outcome measures: change in ESR or DAS28-ESR between baseline and 6 months following methotrexate-based DMARD therapy (statistical analysis by Spearman’s correlation). The fourth (far right) column shows the relationship between the eQTL SNP and the clinical outcome measure: rs2104616, rs6876611 and rs1341486 are linked to change in ESR (right), whereas 12:8206167:A:G is associated with change in DAS28-ESR (statistical analysis by linear model with four genotype PCA covariates). For each of these SNPs, the corresponding eQTL gene (with the exception of ERAP2) is also associated with the response variable. DAS28, Disease Activity Score in 28 joints; DMARD, disease-modifying antirheumatic drug; ESR, erythrocyte sedimentation rate; eQTL, expression quantitative trait locu; PCA, principal component analysis; PEER, Probabilistic Estimation of Expression Residuals.

Figure 3

Selected trans-eQTL effects from leading cis-eSNPs The leading synovial cis-eSNPs were used in a synovial tissue trans-eQTL analysis. Several genes associated with RA disease shared common significant variants. Circos plots, left, show the links between the source SNP and its effect on different trans-eQTL genes across the genome. Boxplots indicate associations between genotype and gene expression (statistical analysis by linear model with four genotype PCA and 4 PEER factors as covariates using matrixEQTL). Scatter plots show correlation between the cis-eQTL gene and trans-eQTL gene (statistical analysis by Spearman’s correlation). (A) SNPs at the ISG15 locus correlated both with baseline ISG15 expression as well as downstream interferon stimulated genes GAS1, IFI44 and IFIT1 and synovial B cell infiltration measured by CD20 histology. (B) SNPs at MZB1 (marginal zone B and B1 cell specific protein) correlated with synovial immunoglobulin gene synthesis and the B-cell related gene CD79A. (C) SNPs at the CDH11 (cadherin 11) locus demonstrated a trans-eQTL correlating with collagen genes and matrix metalloprotease-13 (MMP13). CDH11 trans-eQTL SNPs also correlated with clinical response as demonstrated by change in ESR at 6 months following methotrexate-based DMARD therapy. DMARD, disease-modifying antirheumatic drug; ESR, erythrocyte sedimentation rate; eQTL, expression quantitative trait locu; PCA, principal component analysis; PEER, Probabilistic Estimation of Expression Residuals.

Figure 4

Differences and interactions between synovium and blood eQTL (A) Venn diagram of genome-wide significant genes in synovium and blood with those of known autoimmune disease association in the GWAS catalogue labelled. There are 232 genes which are common to both synovial and blood eQTL. (B) Beta-beta plots highlighting gene-SNP pairs whose eQTLs show an opposite effect in synovium compared with blood. Pairs are coloured according to significance in each tissue and outlined if tissue interaction is significant. (C) Manhattan plot showing the p values for the linear synovial-blood tissue interaction. Significant genes (FDR≤0.01) highlight differences in eQTL between tissues. (D–E) Examples of eQTL effects differing between tissues shown by eQTL locus plots and violin plots for expression in synovium and blood, p values calculated by linear model: NT5E significant tissue interaction; and HIP1 significant in blood only. eQTL, expression quantitative trait locu; FDR, false discovery rate; SNP, single-nucleotide polymorphism.

Figure 5

HLA-DPB2 eQTL is highly associated with disease severity (A–D) Genetic association calculated using linear models via PLINK in HLA region with RA phenotypes (ESR, anti-CCP status, EULAR response and lymphoid subjects). Variants in yellow are imputed HLA residues, those in red are amino acids and those in grey SNPs. Blue markers are the seropositivity markers found by Raychaudhuri et al . (E) Locus plot for association HLA-DPB2 eQTL in synovium using MatrixEQTL with statistical analysis by linear regression with four genotype PCA and four PEER factors as covariates. P values were calculated by linear model t-statistic. Significance plotted on the y axis and ESR significant variants are overlaid with colour-coding. Variants significant as both HLA-DBP2 eSNPS and ESR Candidate Gene Study (CGS) are outlined in black (FDR_eqtl≤0.01, Q_GWAS≤0.01, n=234). (F) Genotype vs HLA-DPB2 expression in synovium for the most significant eQTL variant also significant for ESR (rs3128921). (G–K) Strong correlation (Spearman) between synovial HLA-DPB2 expression with clinical parameters baseline ESR level, baseline disease activity measured by DAS28-ESR, global Visual Analogue Score (VAS) and clinical response to methotrexate-based DMARD treatment assessed by change in ESR between baseline and 6 months. (M–P) P values calculated via linear model with PCA eigenvectors show association between top SNP rs3128921 with the same disease activity and response clinical variables. DAS28, Disease Activity Score in 28 joints; DMARD, disease-modifying antirheumatic drug; eQTL, expression quantitative trait locu; ESR, erythrocyte sedimentation rate; FDR, false discovery rate; GWAS, Genome-Wide Association Studies.