###### Key messages

What is already known about this subject?
=========================================

-   To date, no strong associations of individual genetic loci with response to tumour necrosis factor inhibitors (TNFi) in rheumatoid arthritis (RA) have been identified, despite recent large efforts based on conventional genome-wide association studies and a crowdsourcing initiative.

What does this study add?
=========================

-   We introduced a new methodological approach for localising genetic effects by using genotypic risk scores based on known genetic loci for related traits and likely biomarkers.

-   We identified two genetic loci strongly associated with TNFi response in RA and demonstrated that the genetic determinants of TNFi response are different to the known susceptibility loci for RA.

How might this impact on clinical practice or future developments?
==================================================================

-   Measurements of expression of *CD40* and CD39 and their corresponding pathways could be relevant for targeting drug therapy in RA.

-   Our new methodological approach could be useful for localising genetic effects in traits for which assembling large sample sizes is not feasible, such as drug response.

Introduction {#s1}
============

Biologic therapies have transformed the outlook for rheumatoid arthritis (RA). However, for the most commonly used class of agent, tumour necrosis factor inhibitors (TNFi), there is substantial variability in response to treatment among patients with RA.[@R1] This has spurred efforts to discover predictors of response and more generally to understand how to subtype this heterogeneous disease to predict which therapies will work.[@R2]

Genome-wide association studies (GWAS) of response to TNFi have shown that common single nucleotide polymorphisms (SNPs) explain an estimated 40% and 50% of the variance of change in swollen joint counts (SJC) and erythrocyte sedimentation rate (ESR), respectively; however, no strong associations with individual SNPs have been detected.[@R4] Thus, as with many complex phenotypes, the genetic architecture of response to TNFi is likely to be polygenic with many small genetic effects.[@R5] In this situation, the sample size required to learn a predictive model is very large---up to 10 cases per variable[@R6]---and it may not be feasible to assemble such large sample sizes for studying response to a single drug or drug class.

It has been suggested that improving prediction of complex clinical outcomes may be possible by incorporating information about the genetics of relevant traits in the prediction model.[@R7] One such approach is to use publicly available summary GWAS results of relevant traits to compute genotypic scores, which can then be used as variables ('features') from which to build predictive models. By harnessing the genetic profiles of intermediate traits, these scores aggregate the effects of individual SNPs into larger regional or whole-genome effects. Relevant traits can include diseases, biomarkers and gene transcription levels. For polygenic traits such as RA, for which multiple genetic susceptibility loci have been identified, we can construct locus-specific scores allowing us to examine the extent to which drug response is related to genetic heterogeneity of the disease.

In the current study, we incorporated available genetic information on susceptibility to RA,[@R9] immune cell traits from a publicly available bioresource[@R10] and expression or methylation of genes implicated in response to TNFi treatment in RA.[@R11] The genotypic scores associated with these intermediate traits were then tested for association with response to TNFi; by reducing the number of hypotheses being explored, the thresholds for claiming statistical significance are relaxed, which could help identify useful predictors.

Materials and methods {#s2}
=====================

Cohorts {#s2-1}
-------

We used a sample of 2938 individuals of European ancestry for whom complete clinical and GWAS data were available. This sample comprised individuals from a pre-existing international collaboration formed to study the genetics of response to TNFi agents[@R12] and individuals recruited to the Biologics in Rheumatoid Arthritis Genetics and Genomics Study Syndicate (BRAGGSS) after 2013.[@R4]

[Table 1](#T1){ref-type="table"} shows sample sizes, phenotypes and clinical variables for each of the data collections used in this study. All participants provided informed consent, and institutional review board/ethics approvals were in place as described in Cui *et al* [@R12] and Massey *et al*.[@R4]

###### 

Sample information per cohort.

                                                        BRAGGSS       DREAM         EIRA          ReAct         WTCCC         Other\*       Total
  ----------------------------------------------------- ------------- ------------- ------------- ------------- ------------- ------------- -------------
  Sample size                                           954           764           283           258           556           123           2938
  Sex, female %                                         76            68            74            77            77            82            74
  Concomitant non-biologic DMARD %                      85            74            74            50            73            97            76
  *TNF inihibitor*                                                                                                                          
  Adalimumab                                            416           441           47            258           64            29            1255
  Certolizumab                                          34            0             0             0             0             0             34
  Etanercept                                            293           66            97            0             246           19            721
  Golimumab                                             17            0             0             0             0             0             17
  Infliximab                                            194           138           139           0             246           75            792
  *Baseline disease activity, mean (* *SD* *)*                                                                                              
  DAS28-ESR4                                            6.3 (1.0)     5.4 (1.2)     5.3 (1.2)     5.8 (1.0)     6.7 (0.9)     5.6 (1.1)     6.0 (1.2)
  ESR                                                   36.2 (26.2)   27.7 (21.3)   32.4 (22.9)   31.2 (21.4)   45.0 (28.8)   30.2 (21.5)   34.6 (25.4)
  SJC                                                   10.1 (6.1)    10.3 (5.5)    9.2 (6.0)     9.9 (5.1)     11.8 (6.4)    10.5 (6.3)    10.4 (6.0)
  TJC                                                   15.4 (7.3)    10.0 (7.4)    8.3 (6.0)     13.1 (6.5)    17.0 (7.3)    10.7 (6.3)    13.2 (7.8)
  GHVAS                                                 71.0 (19.2)   62.3 (22.1)   56.1 (23.3)   59.7 (20.8)   72.6 (18.1)   59.3 (23.9)   66.1 (21.4)
  *6* *-* *month* *disease activity, mean (* *SD* *)*                                                                                       
  DAS28-ESR4                                            3.7 (1.6)     3.6 (1.3)     3.5 (1.4)     3.7 (1.4)     4.2 (1.5)     3.9 (1.6)     3.8 (1.5)
  ESR                                                   22.8 (22.0)   18.1 (16.9)   20.0 (17.3)   18.6 (17.0)   27.6 (24.9)   22.2 (19.7)   21.8 (20.7)
  SJC                                                   3.0 (4.0)     4.8 (4.4)     3.4 (3.8)     3.6 (3.6)     4.0 (4.7)     5.0 (4.9)     3.8 (4.3)
  TJC                                                   5.0 (6.3)     3.7 (4.6)     3.6 (4.9)     4.7 (5.5)     6.3 (6.6)     5.1 (6.1)     4.7 (5.8)
  GHVAS                                                 37.2 (25.2)   34.5 (21.7)   34.4 (25.2)   31.8 (25.9)   37.1 (25.0)   31.8 (26.6)   35.5 (24.5)

\* 'Other' displays aggregate sample characteristics for collections with sample size \<100

BRAGGSS, Biologics in Rheumatoid Arthritis Genetics and Genomics Study Syndicate; DAS, Disease Activity Score; DMARD, disease-modifying antirheumatic drug; DREAM, Dutch Rheumatoid Arthritis Monitoring Registry; EIRA, Swedish Epidemiological Investigation of Rheumatoid Arthritis; ESR, erythrocyte sedimentation rate; GHVAS, global health assessment rated on a visual analogue scale; ReAct, French Research in Active Rheumatoid Arthritis; SJC, swollen joint count; TJC, tender joint count; TNF, tumour necrosis factor; WTCCC, Wellcome Trust Case Control Consortium.

Definition of response to TNFi treatment {#s2-2}
----------------------------------------

In RA, response to treatment is quantified by change in the Disease Activity Score (DAS), which depends on four measurements: ESR, SJC, tender joint count (TJC) and patient global health assessment rated on a visual analogue scale (GHVAS). Previous work has shown that only the SJC and ESR measurements have evidence of non-zero heritability[@R4] and correlate significantly with synovitis quantified by ultrasound or MRI.[@R13] Since TNFi were developed to control synovitis, we used the two objective components of the DAS (ESR and SJC) as primary outcomes for evaluating genetic effects and the two subjective components (TJC and GHVAS) and the composite score (DAS28-ESR4) as secondary outcomes. For each outcome, a baseline measurement was taken before initiation of TNFi treatment, and a follow-up measurement was taken between 3 and 6 months after initiation of TNFi treatment. The measurements for each component were transformed in accordance to the DAS28-ESR4 formula (see [online supplementary methods](#SP1){ref-type="supplementary-material"}). Response was modelled as the difference between the baseline and the follow-up measurement.

10.1136/annrheumdis-2018-214877.supp1

Genotypic risk scores {#s2-3}
---------------------

We used the GENOSCORES platform (<https://pm2.phs.ed.ac.uk/genoscores/>) to compute genotypic risk scores for the intermediate traits. GENOSCORES is a database of published SNP to trait associations from a large number of well-powered GWAS, including GWAS of disease traits, biomarkers, gene expression and methylation. The database is accompanied by a software package implemented in R that can be used to compute genotypic risk scores and run downstream statistical analyses in cohorts with SNP data.

We queried the GENOSCORES database for genetic associations with RA risk,[@R9] 149 heritable immune cell traits reported by Roederer *et al*,[@R10] and whole-blood expression and methylation for 93 genes reported in a recent meta-analysis[@R11] as differentially expressed before treatment between responder and non-responder patients with RA treated with TNFi. GWAS summary statistics reported by Westra *et al* [@R15] and Gusev *et al* [@R16] were used for expression quantitative trait loci (eQTLs). GWAS summary statistics reported by Gaunt *et al* [@R17] were used for methylation quantitative trait loci (mQTLs).

GWAS summary statistics for each intermediate trait were filtered at p value\<$10^{- 5}$. SNPs were then split into trait-associated regions, with regions defined as genomic loci containing at least one SNP with p value\<$10^{- 7}$. Only 19 of the immune cell traits had a corresponding trait-associated region. SNPs not assigned to a region were discarded. For each trait-associated region, a genotypic score was computed as a sum of SNP genotypes, $\mathbf{g}$, weighted by the effect size estimates, $\mathbf{\beta}$ (log OR for binary traits, regression slope for quantitative traits) and adjusted for linkage disequilibrium. The regional score, $s_{itr}$, for an individual $i$, trait $t$ and region $r$ was computed as: $s_{itr} = \mathbf{\beta}_{tr}^{T}R_{r}^{- 1}\mathbf{g}_{ir}$, where $R_{r}$ denotes the SNP--SNP correlation matrix in genomic region $r$.

Additional details about the GENOSCORES platform, the score computation and the specific regional scores used in this study are given in online supplementary materials (see [online supplementary methods](#SP1){ref-type="supplementary-material"}, [online supplementary tables S1--S5](#SP2 SP3 SP4 SP5 SP6){ref-type="supplementary-material"}, [online supplementary figure S1](#SP1){ref-type="supplementary-material"}).

Predictive modelling {#s2-4}
--------------------

To evaluate genetic prediction of response to TNFi, we compared a model with clinical covariates only to a model with clinical covariates and genotypic scores for each type of intermediate trait. To avoid numerical instabilities, we removed highly correlated scores prior to fitting a model (see [online supplementary methods](#SP1){ref-type="supplementary-material"}). The number of filtered regional scores for each type of intermediate trait is shown in [table 2](#T2){ref-type="table"}.

###### 

Prediction of response to TNFi using penalised regional genotypic scores for different types of intermediate traits

  -----------------------------------------------------------------------------------------------------------
  Intermediate trait type   No of regional scores   No of filtered scores   Prediction of\   Prediction of\
                                                                            ∆SJC (%)         ∆ESR (%)
  ------------------------- ----------------------- ----------------------- ---------------- ----------------
  Rheumatoid arthritis      37                      37                      5.3 (0.26)       −1.6 (0)

  Immune cell traits        508                     470                     −0.7 (0)         2.9 (0.17)

  eQTLs                     94                      87                      3.4 (0.16)       2.9 (0.17)

  eQTLs and mQTLs           268                     228                     2.5 (0.11)       1.6 (0.09)
  -----------------------------------------------------------------------------------------------------------

Prediction performance is quantified by the difference in test log-likelihood (in nats) between a model with clinical covariates and genotypic scores and a model with clinical covariates only and by the per cent of phenotypic variance explained (in parenthesis). Results from 10-fold cross-validation.

eQTL, expression quantitative trait loci; mQTL, methylation quantitative trait loci.ESR, erythrocyte sedimentation rate; SJC, swollen joint count; TNFi, tumour necrosis factor inhibitors;

We expected that only a subset of genotypic scores would be relevant for prediction of response to TNFi, and thus used a hierarchical shrinkage prior for the score coefficients. We implemented the prediction models in STAN[@R18] using a horseshoe prior distribution and performed inference with Markov chain Monte Carlo sampling.[@R19] To rank the importance of genotypic scores in a model, we applied projection predictive variable selection, an approach that projects posterior draws from the high-dimensional model to lower dimensional subspaces.[@R21]

We used a statistical model with the following clinical covariates: measurements for the four DAS components before initiation of TNFi treatment, gender, whether the patient was concomitantly treated with any non-biologic disease-modifying antirheumatic drugs (DMARDs), cohort (which is also a proxy for country), genotyping array and the 10 first principal components computed from the genotypic data of the full data set. We used individuals with complete measurements in the statistical analyses of each TNFi response outcome (see [online supplementary table S6](#SP1){ref-type="supplementary-material"}).

### Evaluation of prediction {#s2-4-1}

We used two measures to quantify improvement in prediction: the difference in log-likelihood between a model with clinical covariates and genotypic scores and a model with clinical covariates only (measured in natural log units (nats)); and the per cent of residual variance explained by the genotypic scores. Both measures were computed on the testing data from 10-fold cross-validation on the full data set.

For readers who prefer a frequentist interpretation, the asymptotic equivalence of model choice by cross-validation and Akaike's information criterion (AIC)[@R22] implies that a p value of 0.01 for comparison of nested models is equivalent to a difference in test log-likelihood of 2.3 nats (likelihood ratio of 10) for models differing by one extra parameter. The large sample size of this study means that small robust increments in predictive performance can be detected.

Univariate associations {#s2-5}
-----------------------

For models with a test log-likelihood difference of at least two nats, we further examined the univariate associations between genotypic scores and the TNFi response outcomes. We used the full data set to test the univariate association and included the same clinical covariates as in the multivariate prediction, in addition to a score.

For genotypic scores significantly associated with TNFi response at the Bonferroni-corrected p value threshold, we compared the estimated effects among groups receiving different TNFi agents. We considered etanercept, adalimumab and infliximab, where a large enough sample size was available. Additionally, we tested if the associations held when we adjusted for additional covariates: anticitrullinated protein antibodies (ACPA) status and smoking status. These covariates have been reported to influence TNFi response but were only available for a third of the samples, and thus were not included in the full models.

A diagram of the statistical analysis pipeline is given in [online supplementary figure S2](#SP1){ref-type="supplementary-material"}.

Results {#s3}
=======

In the following sections, we focus on the results for the primary TNFi response outcomes. The results for secondary outcomes are discussed in the online supplementary results (see [online supplementary table S7](#SP1){ref-type="supplementary-material"}).

RA genotypic scores {#s3-1}
-------------------

Prediction of response to TNFi as quantified by $\Delta$SJC improved by including the regional genotypic scores for RA risk in a penalised regression model ([table 2](#T2){ref-type="table"}). The test log-likelihood increased by 5.3 nats, suggesting that some of the genetic drivers for RA are also influencing response to TNFi; however, the absolute improvement in prediction was small, with less than 1% of phenotypic variance being explained by the RA genotypic scores.

The regional score at the *CD40* locus had the highest explanatory power for both response phenotypes ([figure 1](#F1){ref-type="fig"}). The direction of the effect was consistent for response as measured by both phenotypes, with higher RA load at the *CD40* locus being associated with better TNFi response. The univariate association of the RA score at the *CD40* locus with $\Delta$SJC passed the p value threshold corrected for the number of RA scores and two response phenotypes ([table 3](#T3){ref-type="table"}).

![Contribution of top 10 RA regional scores to prediction of response to TNFi phenotypes, starting from a model containing only clinical covariates. The curve gradually converges to one with the addition of all remaining scores. ESR, erythrocyte sedimentation rate; RA, rheumatoid arthritis; SJC, swollen joint count; TNFi, tumour necrosis factor inhibitors.](annrheumdis-2018-214877f01){#F1}

###### 

Univariate associations between response phenotypes and regional genotypic scores of interest

  Response phenotype                                              Genetic score           Coefficient   P value
  --------------------------------------------------------------- ----------------------- ------------- ---------
  *Genotypic scores at the CD40 locus*                                                                  
  ∆SJC                                                            RA score at *CD40*      0.07          0.0004
  ∆SJC                                                            *CD40* eQTL             0.06          0.002
  ∆SJC                                                            *CD40* mQTL             −0.05         0.009
  ∆ESR                                                            RA score at *CD40*      0.05          0.01
  ∆ESR                                                            *CD40* eQTL             0.03          0.06
  ∆ESR                                                            *CD40* mQTL             −0.03         0.07
  *Genotypic scores for immune cell traits at the ENTPD1 locus*                                         
  ∆SJC                                                            CD39 on CD4 T           −0.07         5e−05
  ∆SJC                                                            mDC:%32+                −0.02         0.2
  ∆SJC                                                            CD8:%39+                −0.06         0.001
  ∆SJC                                                            CD4:%Treg(39+)          −0.07         0.0001
  ∆SJC                                                            CD4:%Treg(39+73-)       −0.07         0.0001
  ∆SJC                                                            CD4:%Treg(39+73+)       −0.07         0.0002
  ∆SJC                                                            NKeff:%314−158a+        −0.04         0.02
  ∆SJC                                                            CD4 T:%CD39+CD38+PD1−   −0.07         8e−05
  ∆ESR                                                            CD39 on CD4 T           −0.003        0.9
  ∆ESR                                                            mDC:%32+                −0.01         0.5
  ∆ESR                                                            CD8:%39+                −0.008        0.7
  ∆ESR                                                            CD4:%Treg(39+)          −0.004        0.8
  ∆ESR                                                            CD4:%Treg(39+73−)       −0.004        0.8
  ∆ESR                                                            CD4:%Treg(39+73+)       0.006         0.8
  ∆ESR                                                            NKeff:%314−158a+        0.04          0.04
  ∆ESR                                                            CD4 T:%CD39+CD38+PD1−   −0.004        0.8

The coefficients are the effect sizes of the standardised score on the standardised phenotype.

eQTL, expression quantitative trait loci; mTQL, methylation quantitative trait loci.ESR, erythrocyte sedimentation rate; RA, rheumatoid arthritis; SJC, swollen joint count;

The regional score for RA at the *CD40* locus was correlated with a cis-acting eQTL score for *CD40* expression in whole blood (correlation=0.65) and a cis-acting mQTL score for methylation of *CD40* in whole blood (correlation=−0.70). The strongest association between response to TNFi and genotypic scores at the *CD40* locus was with the score for RA ([table 3](#T3){ref-type="table"}). The estimated effect did not change when we stratified by TNFi agent and when we adjusted for ACPA and smoking status (see [online supplementary tables S8 and S9](#SP1){ref-type="supplementary-material"}).

Genotypic scores for immune cell traits {#s3-2}
---------------------------------------

Prediction of response to TNFi as measured by $\Delta$ESR improved by adding penalised genotypic scores for immune cell traits in the model ([table 2](#T2){ref-type="table"}). In univariate analyses, the regional scores for a number of immune cell traits at the *ENTPD1* locus (which codes for CD39) had suggestive associations with $\Delta$SJC ([table 3](#T3){ref-type="table"}). The association of $\Delta$SJC with the genotypic score for 'CD39 on CD4 T cells' at the *ENTPD1* locus passed the p value threshold corrected for the number of immune cell trait scores tested (470) and two response phenotypes. This score was correlated with the cis-acting eQTL score for *ENTPD1* (correlation=0.65). Higher score for 'CD39 on CD4 T cells' at the *ENTPD1* locus was associated with worse TNFi response as quantified by $\Delta$SJC. There was no association between $\Delta$ESR and genotypic scores at the *ENTPD1* locus ([table 3](#T3){ref-type="table"}) and no association between either $\Delta$ESR or $\Delta$SJC and a genotypic score for cell subset frequency of CD73, which is the second ectonucleotidase involved in adenosine production in regulatory T cells (see [online supplementary results](#SP1){ref-type="supplementary-material"}).

Previously, Peres *et al* [@R23] showed that low expression of CD39 on peripheral regulatory T cells was associated with worse response to methotrexate (MTX) in patients with RA. We note that the direction of the effect is reversed, but this is not unlikely since we considered response to TNFi. To further investigate the effect of 'CD39 on CD4 T cells' expression on TNFi response, we performed two analyses stratified by concomitant treatment. Information on whether a patient was receiving a concomitant non-biologic DMARD was available for all samples, while information on whether a patient was specifically receiving concomitant MTX treatment was available for a subset of patients from the BRAGGSS cohort. [Table 4](#T4){ref-type="table"} shows the effect of the genotypic score for expression of 'CD39 on CD4 T cells' on $\Delta$SJC for patients receiving TNFi treatment stratified by concomitant treatment with either any non-biologic DMARD (top) or specifically with MTX (bottom). The effect became stronger in the groups receiving concomitant treatment and attenuated in the group receiving TNFi monotherapy; the CIs among all groups overlapped. Similarly, we did not detect statistically significant differences when we stratified by TNFi agent and when we adjusted for ACPA status (see [online supplementary tables S8 and S9](#SP1){ref-type="supplementary-material"}).

###### 

Univariate association between ∆SJC and genotypic score for the expression of 'CD39 on CD4 T cells' at the *ENTPD1* locus stratified by concomitant treatments

  Patient group                                 Coefficient (SE)   P value   Sample size
  --------------------------------------------- ------------------ --------- -------------
  All samples, adjusted for concomitant DMARD   −0.07 (0.02)       5e−05     2922
  All samples, no adjustment for DMARD          −0.08 (0.02)       2e−05     2922
  Samples not on concomitant DMARD              −0.02 (0.04)       0.5       691
  Samples on concomitant DMARD                  −0.09 (0.02)       1e−05     2231
  Samples on concomitant MTX                    −0.1 (0.03)        0.002     958

The coefficients are the effect sizes of the standardised score on the standardised phenotype. The phenotype is adjusted for covariates: baseline DAS components, gender, cohort, genotyping array, 10 genetic principal components.

DAS, Disease Activity Score; DMARD, disease-modifying antirheumatic drug; MTX, methotrexate; SJC, swollen joint count.

eQTL and mQTL scores {#s3-3}
--------------------

Prediction of response to TNFi as quantified by both $\Delta$SJC and $\Delta$ESR improved by adding eQTL and mQTL scores of implicated genes in a penalised regression model ([table 2](#T2){ref-type="table"}). Of the 93 genes reported in Kim *et al* [@R11] as differentially expressed between responders and non-responders to TNFi in RA, 54 genes had at least one eQTL, 54 genes had at least one mQTL and 36 had both. The test log-likelihood increased by 3.4 and 2.9 nats for $\Delta$SJC and $\Delta$ESR, respectively, by adding the eQTL scores in the model. We did not see a further improvement by adding genotypic scores for mQTLs.

Discussion {#s4}
==========

In the largest international study of TNFi response to date, we have shown how using methods that leverage information from relevant intermediate traits can identify predictors of TNFi response. In a recent crowdsourced effort to use machine learning to construct a predictor of response to TNFi in RA, including SNP genotypes did not improve prediction beyond that obtained with clinical covariates alone.[@R24] Genotypic prediction of psychiatric disorders and related phenotypes has been shown to improve by exploiting genetic correlations among multiple related traits,[@R7] and methods have been extended to incorporate polygenic scores for multiple traits.[@R8]

In the current study, we have combined these approaches and implemented them in a newly developed platform, called GENOSCORES, which contains GWAS data for multiple traits and automates construction of genotypic scores. For polygenic traits with multiple trait-associated loci, such as RA, locus-specific scores can be constructed to examine how genetic heterogeneity of the intermediate trait can influence the trait of interest. Our approach reduces the dimensionality of the prediction task from about 2 million common SNPs to a few hundred or a few thousand genotypic scores, depending on how relevant traits are selected. The score constructions are a type of feature engineering, a task commonly used in machine learning applications.

Understanding the pathogenic mechanisms that initiate and perpetuate RA could give rise to informative biomarkers of prognosis, therapeutic response and toxicity.[@R26] However, in agreement with earlier studies,[@R27] we did not find strong predictors of TNFi response among alleles linked to the development of RA. A strength of the current study is the large sample size which allowed us to detect small robust increments in predictive performance. Our methodological approach was to first establish the predictive value for a set of genetic markers using a multivariate model and then to examine univariate associations between each marker and the outcome. Using this approach, we showed that a model including RA scores led to a small robust improvement in prediction, with the regional score at the *CD40* locus driving the predictive signal.

Higher RA risk at the *CD40* locus, higher *CD40* transcription and lower *CD40* methylation were associated with better TNFi response. CD40 is a transmembrane protein which belongs to the TNF receptor superfamily, critically important in modulating immune-(auto-immune) responses.[@R28] CD40 is expressed by B cells and antigen-presenting cells (APCs), whereas CD40L is induced on CD4+ T cells following T-cell antigen receptor (TCR) with major histocompatability complex (MHC) molecule interaction. Engagement of the CD40--CD40L axis leads to B cell activation, proliferation and (auto)-antibody production, while activation of APCs by CD40L on CD4+ T cells induces upregulation of CD80, CD86, MHC class I and MHC class II, as well as secretion of proinflammatory cytokines such as interleukin (IL)-12, IL-23 and TNF-$\alpha$.[@R29]

The risk allele associated with RA is associated with elevated *CD40* expression in whole blood.[@R31] As high *CD40* expression is associated with elevated TNF-$\alpha$ production and *CD40* and *CD40L* transcripts are increased in the disease tissue in both early and established disease,[@R32] it is not surprising that patients with the *CD40* risk allele respond better to TNFi therapies.

Overall, if there are genetic loci that predict TNFi response, these are mostly different to the known RA risk loci. However, we note that patients who receive TNFi therapy are likely to have more severe disease and to have failed on other treatments. It is therefore possible that our study sample has been selected with respect to genetic load for RA, thus limiting the heterogeneity in genetic RA risk profiles compared with a sample of newly diagnosed cases.

The genotypic score for the expression of the ectonucleotidase CD39 on CD4 T cells was inversely associated with TNFi response. The SNPs contributing to this score are in the *ENTPD1* gene which encodes CD39. In stratified analyses, the inverse association with response was stronger in the groups receiving TNFi concomitantly with MTX or another non-biologic DMARD compared with the group receiving TNFi monotherapy and was stronger in the group receiving infliximab compared with the groups receiving adalimumab or etanercept. The CIs of the estimated effects among all groups overlapped. This effect on response to TNFi agents is in the opposite direction to the association reported between low expression of CD39 on regulatory T cells and resistance to MTX in RA.[@R23]

Interpreting the association between drug response and the CD39 trait is difficult both epidemiologically and mechanistically. The reasons for this being that in the UK and most European countries TNFi are usually prescribed only after patients have had a poor response to MTX and, unless of intolerance, always in combination with MTX. Therefore, these patients are likely to represent a selected group; this, together with the almost universal use of combination MTX/TNFi therapy, hinders progress in dissecting the potential mechanisms responsible for the divergence. Nonetheless, as RA is a highly heterogeneous disease, it is plausible to speculate that different cellular and molecular networks may be involved in driving diverse immune/inflammatory responses in different patients to different drugs. For example, as mentioned above, the poor response to MTX has been associated with low CD39 expression by regulatory T cells, while increased CD39 expression has been reported to be important in the expansion of Th17 cells driven by IL-6 and TGF-β via Stat3 and Gfi-1 transcription factors.[@R33] In turn, the expansion of Th17 cells has been reported to be associated with incomplete response to TNFi.[@R34] It remains to be established whether measurements of CD39/*ENTPD1* expression or genotype may be useful in the choice of MTX, TNFi agent or concomitant treatment as first-line therapy for patients who need a DMARD.

Using previously reported associations between transcripts and TNFi response to select relevant genes, we have shown evidence that eQTL scores for these genes contain information that predicts TNFi response, even though the proportion of variance explained was low and no single genes associated with response could be identified.

Improved genomic prediction of treatment response requires measuring response more precisely to capture the molecular in addition to the clinical phenotype. The detected associations between genotypic scores and TNFi response were not the same for the two measures of response---change in SJC and change in ESR---suggesting that the two measures reflect different aspects of disease activity affected by TNFi. To derive refined measures of drug response, large data sets with multiple inflammatory biomarkers, joint imaging and clinical variables before and after treatment are needed.

10.1136/annrheumdis-2018-214877.supp2

10.1136/annrheumdis-2018-214877.supp3

10.1136/annrheumdis-2018-214877.supp4

10.1136/annrheumdis-2018-214877.supp5

10.1136/annrheumdis-2018-214877.supp6

10.1136/annrheumdis-2018-214877.supp7

We thank all the patients who have contributed to this research and clinical staff who supported patient recruitment. We thank the Medical Research Council and Arthritis Research UK (ARUK) for their joint funding of the Maximising Therapeutic Utility in Rheumatoid Arthritis (MATURA) study. We thank MATURA consortium members for their input and feedback during the design of this project, and we would like to give particular mention to Felix Agakov, Jennifer H Barrett and Heather J Cordell. The list of MATURA collaborators and affiliations is available in the online supplementary information (table S12). We thank Peter K Gregersen for providing access to clinical and genetic data from the ABCoN cohort. We thank Stephanie Ling for collating ACPA data from the BRAGGSS cohort into a usable format. We thank the anonymous reviewers for their constructive comments and suggestions.

**Handling editor:** Josef S Smolen

**Correction notice:** This article has been corrected since it first published. Table 1 data aligmment has been corrected.

**Collaborators:** MATURA Consortium.

**Contributors:** Substantial contributions to conception or design of the study: AS, PM, AB, CP. Substantial contributions to drafting the manuscript: AS, PM. Substantial contributions to data acquisition: JC, MJHC, KI, HY, SS, LP, SLB, RPK, YO, PLCMvR, GW, IEvdH-B, NdV, PPT, KO, HC, H-JG, TWJH, LAC, SR, MW, AGW, XM, JDI, AWM, AB. Substantial contributions to data analysis or interpretation: AS, PM, MC, DP, NN, CP. All authors contributed to revising the manuscript critically for important intellectual content and approved the final manuscript. The funding agencies had no part in writing or reviewing the manuscript.

**Funding:** This study was supported by Medical Research Council (MRC) MR/K015346/1 MATURA study. ARUK 20670 MATURA study.

**Competing interests:** None declared.

**Patient and public involvement statement:** This study was conducted as part of Work Stream 2 of the MATURA Consortium (<http://www.matura.whri.qmul.ac.uk/what_is_matura.php>). A patient advisory group was established in 2014 when the MATURA project commenced. The group meets regularly to:

Ensure MATURA strategy is maintaining relevance, accountability and direction by embedding patients and members of the public within the decision making processes.

Determine what level of confidence in tests, and what type of tests, would be acceptable to patients for treatment decisions.

Maximise patient recruitment to research studies by increasing awareness through patient groups.

Readily obtain patients' perspective on grant applications related to stratified medicines for RA.

Facilitate the dissemination of the results from MATURA research, for instance by producing lay summaries of papers in conjunction with the researchers (<http://www.matura.whri.qmul.ac.uk/news.php>).

**Patient consent for publication:** Not required.

**Ethics approval:** This study used anonymised data for human subjects from an international collaboration of 13 studies originally published in PLOS Genetics 2013;9:e1003394. All participants provided informed consent and institutional review board and ethics approvals were in place for each of the studies and described in the original publication.

**Provenance and peer review:** Not commissioned; externally peer reviewed.

**Data sharing statement:** Data are available upon reasonable request.
