Open Access

A longitudinal genome-wide association study of anti-tumor necrosis factor response among Japanese patients with rheumatoid arthritis

  • Kyoko Honne1,
  • Ingileif Hallgrímsdóttir2,
  • Chunsen Wu3, 4,
  • Ronnie Sebro5,
  • Nicholas P. Jewell6,
  • Takeo Sakurai7,
  • Masahiro Iwamoto1,
  • Seiji Minota1 and
  • Damini Jawaheer8Email author
Arthritis Research & Therapy201618:12

https://doi.org/10.1186/s13075-016-0920-6

Received: 9 June 2015

Accepted: 4 January 2016

Published: 18 January 2016

Abstract

Background

Studies of Caucasian patients with rheumatoid arthritis (RA) to identify genetic biomarkers of anti-tumor necrosis factor (TNF) response have used response at a single time point as the phenotype with which single nucleotide polymorphism (SNP) associations have been tested. The findings have been inconsistent across studies. Among Japanese patients, only a few SNPs have been investigated. We report here the first genome-wide association study (GWAS) to identify genetic biomarkers of anti-TNF response among Japanese RA patients, using response at 2 time-points for a more reliable clinical phenotype over time.

Methods

Disease Activity Scores based on 28 joint counts (DAS28) were assessed at baseline (before initial therapy), and after 3 and 6 months in 487 Japanese RA patients starting anti-TNF therapy for the first time or switching to a new anti-TNF agent. A genome-wide panel of SNPs was genotyped and additional SNPs were imputed. Using change in DAS28 scores from baseline at both 3 (ΔDAS-3) and 6 months (ΔDAS-6) as the response phenotype, a longitudinal genome-wide association analysis was conducted using generalized estimating equations (GEE) models, adjusting for baseline DAS28, treatment duration, type of anti-TNF agent and concomitant methotrexate. Cross-sectional analyses were performed using multivariate linear regression models, with response from a single time point (ΔDAS-3 or ΔDAS-6) as phenotype; all other variables were the same as in the GEE models.

Results

In the GEE models, borderline significant association was observed at 3 chromosomal regions (6q15: rs284515, p = 6.6x10−7; 6q27: rs75908454, p = 6.3x10−7 and 10q25.3: rs1679568, p = 8.1x10−7), extending to numerous SNPs in linkage disequilibrium (LD) across each region. Potential candidate genes in these regions include MAP3K7, BACH2 (6q15), GFRA1 (10q25.3), and WDR27 (6q27). The association at GFRA1 replicates a previous finding from a Caucasian dataset. In the cross-sectional analyses, ΔDAS-6 was significantly associated with the 6q15 locus (rs284511, p = 2.5x10−8). No other significant or borderline significant associations were identified.

Conclusion

Three genomic regions demonstrated significant or borderline significant associations with anti-TNF response in our dataset of Japanese RA patients, including a locus previously associated among Caucasians. Using repeated measures of response as phenotype enhanced the power to detect these associations.

Keywords

Rheumatoid arthritis Anti-TNF GWAS Japanese Pharmacogenomics Pharmacogenetics MAP3K7 BACH2 GFRA1 WDR27

Background

Anti-tumor necrosis factor (TNF) agents have been particularly effective in controlling disease activity and joint erosions in rheumatoid arthritis (RA) [13]. Nonetheless, a considerable proportion of patients (30–40 %) demonstrate only partial response or no response to these therapies. This complex phenotype of anti-TNF responsiveness has a genetic component, as demonstrated by heritability estimates of 0.59 or more for different clinical measures of response [4]. Several studies conducted mainly in Caucasian populations [523] have attempted to identify genetic biomarkers that can reliably predict response to anti-TNF agents in patients with RA. However, there have been no reports of genome-wide association studies (GWAS) for biomarkers of anti-TNF response among Japanese or other East Asian RA patients. Only single nucleotide polymorphisms (SNPs) in the TRAF1 [24] and CD84 genes [18] have been examined for association in two small samples of 101 and 151 Japanese patients, respectively. It has been demonstrated for several RA risk loci that there are both similarities and ethnic differences in disease associations between Caucasian and Japanese populations [2529]. Similarly, although efficacy of anti-TNF therapies appears to be similar in these populations [30, 31], there may be similarities and differences in genetic predictors of anti-TNF response.

Previous studies investigating genetic predictors of anti-TNF response focused on a limited number of candidate genes [515, 2022, 32], and few GWAS have been performed [9, 1619, 23]. The findings from these different studies have been largely inconsistent. Only a couple of loci (PTPRC and PDE3A-SLCO1C1) have been associated in more than one study [5, 11, 33], but these associations were not consistently replicated in other studies [8, 14, 20]. Findings at other loci have not been independently replicated [12, 22]. All of these studies used a single assessment of anti-TNF response.

In the present study, we have performed a GWAS of anti-TNF response in a population of Japanese RA patients. In contrast to previous studies, we have used repeated measures of response at 3 and 6 months after initiation of anti-TNF therapy in order to have a more reliable assessment of clinical phenotype over time, and have used longitudinal statistical models to accommodate the repeated measures of the outcome.

Methods

Patients

A total of 487 RA patients were recruited from the Jichi Medical University and from Inoue Hospital within the Gunma prefecture in Japan, and prospectively followed for 6 months. All included patients: (1) were of Japanese descent; (2) satisfied the 1987 American College of Rheumatology (ACR) criteria for RA [34]; (3) were not in remission (Disease activity score in 28 joints (DAS28) <2.6) at baseline; and (4) were starting treatment with the anti-TNF agents etanercept (ETN), infliximab (INF) or adalimumab (ADA) for the first time or were switching to a new anti-TNF drug due to inefficacy or adverse events.

Data and sample collection

Study baseline

Detailed clinical and demographic data, including age, sex, disease duration, smoking history, and concurrent treatment with methotrexate and prednisolone, were collected on each patient prior to initiation of anti-TNF therapy (study baseline). A serum sample was collected, and used to measure baseline titers of anti-citrullinated protein antibodies (ACPA), rheumatoid factor (RF) using the RA particle agglutination test (RAPA), and C-reactive protein (CRP) levels. Patients were considered to be seropositive for ACPA and/or RF auto-antibodies if ACPA and RAPA titers were ≥4.5 U/ml and ≥40 units, respectively. Whole blood samples (8.5 ml) were collected from each patient and stored at −80 °C for subsequent DNA extraction.

Follow up

All patients were evaluated at 3 and 6 months after initiation of anti-TNF therapy. Tender joint counts (TJC28) and swollen joint counts (SJC28) for 28 joints, patient global scores and CRP levels were assessed at each of these time points. The study protocol was approved by the ethics committee of the Jichi Medical University and by the Institutional Review Board of the Children’s Hospital Oakland Research Institute. Signed informed consent was obtained from each study participant.

Assessment of response to anti-TNF therapy

RA disease activity was assessed at baseline, 3 months and 6 months by the DAS28 calculated on three variables, including 28-joint counts and CRP, i.e., DAS28CRP3 [35, 36], as follows:
$$ \mathrm{D}\mathrm{A}\mathrm{S}28\mathrm{C}\mathrm{R}\mathrm{P}3=\left(0.56\ast \sqrt{TJC28}\right)+\left(0.28\ast \sqrt{SJC28}\right)+\left(0.36\ast ln\left(CRP+1\right)\ast 1.10\right)+1.15 $$

We will refer to the DAS28CRP3 as DAS28. The change in DAS28 scores from baseline at 3 months and 6 months were calculated as ΔDAS-3 and ΔDAS-6, respectively.

Genotyping and data cleaning

DNA was extracted from frozen whole blood samples using a standard protocol and a Gentra Autopure system (Gentra Systems, Minneapolis, MN, USA). All samples were genotyped for a total of 1,133,484 SNPs on an Illumina BeadLab1000 platform using the HumanOmni1-Quad BeadChip and the Infinium HD assay (Illumina, San Diego, CA, USA). Genotypes were called using a score threshold of 0.15 in the Illumina BeadStudio software. Data cleaning was performed using the PLINK software [37] (http://pngu.mgh.harvard.edu/~purcell/plink/) SNPs with genotyping rates ≤98 %, minor allele frequencies (MAF) ≤1 %, or not in Hardy Weinberg equilibrium (HWE) (p ≤0.0001) were excluded from further analyses due to possible genotyping error. To test for population stratification, principal component analysis (PCA) was performed using the EIGENSTRAT software [38, 39].

Imputation of genotypes

The cleaned genotypes were phased using the ShapeIT (v2) software [40] and imputation of genotypes was performed using the Impute2 software [41]. All available multi-population haplotypes from the 1000 Genomes haplotypes Phase I integrated variant set (June 2014 release) were used as reference panels both for phasing and imputation, as recommended [42]. The probability distribution of three possible genotypes generated by Impute2 at each imputed SNP was converted to genotypes using the GTOOL software (v 0.7.5) (http://www.well.ox.ac.uk/~cfreeman/software/gwas/gtool.html) and a stringent probability threshold of 0.9 was applied. Imputed SNPs with genotyping rates ≤98 % or MAF <5 % were excluded from subsequent analyses.

Statistical analyses

Variables influencing response to anti-TNF therapy

To identify potential confounder variables that influence patient response to anti-TNF therapy over time, a longitudinal analysis was performed using GEE models to accommodate response at two time-points for each patient and to adjust for within-patient correlation [43]. Repeated measures of the change in DAS28 at 3 and 6 months (i.e., ΔDAS-3 and ΔDAS-6) were used as the outcome variable in the models, and the explanatory variables included baseline DAS28, duration of anti-TNF therapy, age at baseline, RA duration, sex, concurrent methotrexate use (yes/no), concurrent prednisolone use (yes/no), type of anti-TNF agent (ETN, INF or ADA), RAPA (yes/no), smoking status (never/ever) and ACPA seropositivity (yes/no). Each variable was tested for association with repeated measures of the change in DAS28 in univariate and multivariate models.

Longitudinal genome-wide association analyses

Longitudinal GEE models were used to investigate associations between each SNP on a genome-wide scale and patient response to anti-TNF therapy over time (response at two time points for each patient). For each SNP, repeated measures of the change in DAS28 at 3 and 6 months (i.e., ΔDAS-3 and ΔDAS-6) were used as the outcome variable in the model, with the SNP being the explanatory variable. Covariates included in the model were baseline DAS28, concurrent methotrexate use (yes/no), type of anti-TNF agent (ETN, INF or ADA) and duration of anti-TNF therapy. These covariates were selected on the basis of their association with change in DAS28 from the GEE models described earlier. Principal components were not adjusted for in the main model because the genomic control inflation factor (λGC) was estimated at 1.001. These analyses were repeated after excluding patients who had mild disease activity at baseline.

All analyses were performed using the STATA package (version 13). P values below a non-stringent threshold of 1×10−6 were taken as evidence of borderline significant association.

Cross-sectional genome-wide association analyses at each time point

Cross-sectional genome-wide association (GWA) analyses were performed using response data from a single time point, at 3 months or at 6 months, as the phenotype in two separate multivariate linear regression models. These models (model 1 and model 2) differed from the GEE models only in terms of the response phenotype: ΔDAS-3 was the outcome variable in model 1, and ΔDAS-6 in model 2. All other explanatory variables and covariates were the same as those in the GEE models. The GWA analyses for models 1 and 2 were performed using the PLINK software.

Identifying regions of association

Independent chromosomal regions of association were identified using the LD clumping option within PLINK, based on p values and patterns of LD in the data.

Results

Patient characteristics

A total of 444 patients had genotype data and data on ΔDAS-3, ΔDAS-6, and covariates available for analysis. Among these, ΔDAS-3 was missing for 3 patients and ΔDAS-6 was missing for 22 patients. The clinical, demographic and treatment characteristics of the patients at baseline are summarized in Table 1. The majority of patients were female (84 %), with RA duration of 8.1±8.5 years, had moderate to severe disease activity at baseline (93 %), were naïve to anti-TNF drugs (94 %), were concurrently being treated with methotrexate (80 %) and were positive for ACPA antibodies (89 %). Overall, response to anti-TNF therapy was significantly better at 6 months than at 3 months (ΔDAS-3 (mean ± SD): 1.50±1.01, ΔDAS-6: 1.72±1.12; p = 0.002).
Table 1

Patient characteristics at baseline

Characteristic

All patients (n = 444)

Age, years

56.4 ± 12.7

Sex, female

373 (84.0 %)

RA duration (years)

8.1 ± 8.5

Serum CRP (mg/l)

27.6 ± 27.9

Tender joint count

6 (3–11)

Swollen joint count

6 (3–10)

DAS28CRP3 score

4.6 ± 1.0

Disease activity

 

 Mild (2.6 < DAS28 ≤ 3.2)

31 (7.0 %)

 Moderate (3.2 < DAS28 ≤ 5.1)

282 (63.5 %)

 Severe (DAS28 > 5.1)

131 (29.5 %)

Anti-TNF therapy

 

 Etarnercept

172 (38.7 %)

 Infliximab

242 (54.5 %)

 Adalimumab

30 (6.8 %)

Other medications

 

 Methotrexate

354 (79.7 %)

 Prednisolone

283 (65.1 %)*a

Anti-TNF naïve at baseline

419 (94.4 %)

Switching to a new anti-TNF agent at baseline

18 (4.1 %)

Ever smoker

94 (23.4 %)*b

ACPA positive

330 (88.7 %)*c

RF positive

337 (80.1 %)*d

Results are shown as mean ± SD, median (IQR) or number (%). *Data for these variables were not available for all patients. Sample sizes were as follows a435, b402, c372, d421. RA rheumatoid arthritis, CRP C-reactive protein, Disease Activity Score in 28 joints DAS28, ACPA anti-citrullinated protein antibodies, RF rheumatoid factor

Variables influencing response to anti-TNF therapy

In univariate and multivariate GEE models, only baseline DAS28 (p <0.0005), concurrent methotrexate use (p <0.0005) and duration of anti-TNF therapy (p <0.0005) were associated with change in DAS28 (Table 2).
Table 2

Results from multivariate GEE model to identify predictors of response

Variable

β (95 % CI)

P value

Age at baseline, years

−0.003 (−0.01, 0.004)

0.35

RA duration, years

0.001 (−0.01, 0.008)

0.84

Sex, malea/female

−0.13 (−0.41, 0.16)

0.39

DAS28 at baseline

0.39 (0.30, 0.49)

<0.0005

Duration of anti-TNF therapy, weeks

0.02 (0.01, 0.03)

<0.0005

RAPA seropositivity, yes/noa

−0.07 (−0.39, 0.25)

0.66

ACPA seropositivity, yes/noa

0.10 (−0.28, 0.48)

0.61

Smoking status, nevera/past/current

−0.10 (−0.24, 0.05)

0.18

Concurrent methotrexate use, yes/noa

0.43 (0.19, 0.68)

<0.0005

Concurrent prednisolone use, yes/noa

−0.08 (−0.28, 0.11)

0.40

Type of anti-TNF therapy, ETNa, ADA or INF

−0.08 (−0.18, 0.02)

0.12

β regression coefficient. aReference category for these variables in the generalized estimating equations (GEE) model, RA rheumatoid arthritis, DAS28 Disease Activity Score in 28 joints, RAPA rheumatoid arthritis particle agglutination, ACPA anti-citrullinated protein antibodies

Data quality and population stratification

A total of 4,253,138 autosomal SNPs were available for statistical analyses after imputation and after quality control thresholds were applied. Of these, 738,576 had been genotyped, and the rest had been imputed with high accuracy (imputation accuracy score >0.96).

SNP associations

GEE models using repeated measures of change in DAS28 as the phenotype

Three independent genomic regions showed borderline significant (p <1×10−6) association with repeated measures of change in DAS28 (i.e., both ΔDAS-3 and ΔDAS-6) (Fig. 1); the results for the index SNPs within each region are shown in Table 3. Two of these regions mapped to chromosome 6: one at 6q15 (rs284515: p = 6.6×10−7) approximately 15 Kb downstream from the Mitogen-activated protein kinase kinase kinase 7 (MAP3K7) gene and 572 Kb upstream from the Basic leucine zipper transcription factor 2 (BACH2) gene, and the other at 6q27 (rs75908454: p = 6.3×10−7) within a locus containing the WD repeat-containing protein 27 (WDR27) gene. The third locus (rs1679568: p = 8.1×10−7) mapped to the 3’ untranslated region (UTR) of the Glial cell line-derived neurotrophic factor family receptor alpha 1 (GFRA1) gene at 10q25.3. As shown in Fig. 1, at each locus, numerous SNPs within the region of LD containing the index SNP also showed evidence of association with change in DAS28. Similar results were obtained when the analysis was restricted to patients with moderate or severe disease activity at baseline (n = 413) (see Additional file 1). In addition, the results were again similar when the analysis was restricted only to the subset of patients who were anti-TNF naïve at baseline (n = 419). In both cases though, the significance levels were lower than in the full dataset of 444 patients. Chromosomal regions showing moderate evidence of association (p <1×10−5) are shown in Additional file 2, and the index SNPs within those regions are listed in Additional file 3.
Fig. 1

Regions demonstrating significant or borderline significant association with anti-TNF response. The regional plots show three regions with borderline significant association (p <1×10−6) at 6q15 (top row), 6q27 (middle row) and 10q25.3 (bottom row), based on results from all patients (n = 444) and generalized estimating equations (GEE) models using anti-TNF response at two time points, i.e., 3 and 6 months (left panels). Association results at the same genomic regions are shown for models using anti-TNF response at a single time point only: 6 months (middle panels) and 3 months (right panels). Note: at 6 months, another single nucleotide polymorphism (SNP) in linkage disequilibrium with rs284515 the 6q15 locus (i.e., rs284511) was significantly associated with response at 6 months (change in Disease Activity Score at 6 months, ΔDAS-6); this SNP was not associated at 3 months, and was borderline significant in the GEE model. Chr chromosome

Table 3

Index SNPs from independent chromosomal regions showing borderline significant association with repeated measures of response at 3 and 6 months

SNP

Chr

Position (bp)

Minor allele

MAF

Genotype

N

β (95 % CI)

P value

Gene(s) in region of association

rs284515

6

91,208,828

G

0.17

AG/GG

134

0.35 (0.21, 0.49)

6.6×10−7

MAP3K7, BACH2

     

AG

121

0.33 (0.16, 0.51)

2.1×10−4

     

GG

13

0.80 (0.43, 1.18)

2.6×10−5

rs75908454

6

169,970,623

C

0.10

CT/CC

88

0.42 (0.26, 0.59)

6.3×10−7

WDR27

     

CT

87

0.41 (0.24, 0.59)

2.3×10−6

     

CC

1

-

-

rs1679568

10

117,817,551

A

0.15

AG/AA

124

0.35 (0.21, 0.49)

8.1×10−7

GFRA1

     

AG

114

0.37 (0.20, 0.54)

1.6×10−5

     

AA

10

0.63 (0.24, 1.02)

1.6×10−3

The single nucleotide polymorphisms (SNP) with the smallest p value (index SNP) in each of the associated regions are shown. β (regression coefficient) represents the increase in change in DAS28 associated with the genotype(s) shown compared to the reference genotype (with no minor alleles). Chr chromosome, MAF minor allele frequency, N number of patients with genotype(s) shown

The CD84, PTPRC, PDE3A-SLCO1C1 and MED15 loci reported to be significantly associated with anti-TNF response in RA in previous GWAS showed no evidence of association in our dataset.

Cross-sectional analyses at 3 months and at 6 months

When response at only the 3-month or 6-month follow up (ΔDAS-3 or ΔDAS-6) was used as the phenotype instead of ΔDAS at both time points, SNP rs284511 mapping close to and in LD with rs284515 near the MAP3K7 locus was significantly associated with ΔDAS-6 at the genome-wide significance level (p = 2.5×10−8). No other SNPs showed significant or borderline significant associations with either ΔDAS-6 or ΔDAS-3.

Of interest, as shown in Fig. 1 and Table 4, compared to the cross-sectional analyses which used response phenotype from a single time point, the GEE models provided increased power to detect associations at loci that demonstrated modest association at both time points, i.e., rs284515, rs75908454 and rs1679568, but not rs284511.
Table 4

SNPs showing significant or borderline significant association in the cross-sectional analyses or in the GEE model

     

Model 1

Model 2

GEE model

SNP

Chr

Position (bp)

Minor allele

MAF

Phenotype: ΔDAS-3

Phenotype: ΔDAS-6

Phenotype: ΔDAS-3 and ΔDAS-6

N

β

P

N

β

P

N

β

P

rs284511

6

91,208,542

C

0.33

440

0.16

0.01

421

0.41

2.5×10−8

443

0.28

1.5×10−6

rs284515

6

91,208,828

G

0.17

438

0.35

5.6×10−5

419

0.36

1.1×10−4

443

0.35

6.6×10−7

rs75908454

6

169,970,623

C

0.10

440

0.36

0.0007

421

0.49

5.5×10−5

443

0.42

6.3×10−7

rs1679568

10

117,817,551

A

0.15

434

0.27

0.002

415

0.44

5.3×10−6

437

0.35

8.1×10−7

Association results from cross-sectional (models 1 and 2) and longitudinal (generalized estimating equations (GEE) model) analyses - for single nucleotide polymorphisms (SNPs) within regions showing borderline significant association with the GEE model. MAF minor allele frequency ΔDAS-3 and ΔDAS-6, change in Disease Activity Scores based on 28 joint counts (from baseline) at 3 months and 6 months, Chr chromosome

Discussion

This is the first GWAS investigating genetic biomarkers of response to anti-TNF therapy in Japanese patients with RA, utilizing a longitudinal approach to examine associations between genome-wide SNPs and repeated measures of anti-TNF response at 3 and 6 months. We found borderline significant association (p <1×10−6) at three non-correlated regions within our study population, with the associated SNPs mapping to or close to the following genes: MAP3K7, BACH2 (6q15), WDR27 (6q27) and GFRA1 (10q25.3). Each of these regions harbored numerous SNPs demonstrating evidence of association with change in DAS28. Furthermore, the 6q15 locus was significantly associated with response at 6 months (p = 2.5×10−8), and the association at the GFRA1 locus represents a replication of a previously reported association among Caucasian patients. We therefore considered these regions worthy of being reported so that they may be investigated further in larger datasets.

The MAP3K7 gene encodes transforming growth factor beta-activated kinase 1 (TAK1) which is a key regulator in multiple inflammatory signaling pathways [44, 45], including the p38 MAPK and nuclear factor kappa B signaling pathways. TAK1 deficiency leads to reduced pro-inflammatory cytokine production in cultured RA synoviocytes [46]. It is thus an excellent candidate that may influence the effect of anti-TNF agents, as proposed [47], and is already a candidate therapeutic target to block pro-inflammatory pathways in RA [48, 49]. Transcription factor BACH2, on the other hand, appears to be a key negative regulator of effector T cell differentiation, promoting immune homeostasis [50]. In the mouse, it appears to be a super-enhancer repressing a network of genes critical for T cell function [51]. Of interest, variants at the BACH2 locus have been associated with multiple autoimmune diseases, including RA [5257]. The GFRA1 protein is a member of the Glial cell line-derived neurotrophic factor (GDNF) receptor family and mediates activation of the RET tyrosine kinase receptor. GDNF is produced by astrocytes in response to pro-inflammatory cytokines including TNFα [58] and appears to suppress interleukin-17 (IL-17)-mediated inflammation via the NF-kappa B pathway [59]. The function of the WDR27 protein has not been established.

The association with the GFRA1 gene was previously identified by Plant et al. [16]. SNP rs7070180 mapping to an intron of the GFRA1 gene, was associated with anti-TNF response in a cohort of 566 Caucasian RA patients in the UK (p = 2.24×10−4) and in a meta-analysis including a cohort of 379 additional patients (p = 6.42×10−5). However, this locus was not reported among the major findings of that study as SNP rs7070180 failed to genotype in one of the cohorts. While this SNP was not genotyped or imputed in our data, several other SNPs within the 3’ untranslated region (UTR) of the GFRA1 gene were associated with anti-TNF response among our patients. There are no reports of MAP3K7 being associated with anti-TNF response in RA, although it has been proposed as a good candidate for pathway pharmacogenetics relating to TNF inhibitors [47]. Among the genes of the p38 MAPK network that have been investigated [7, 9, 13], suggestive evidence of an association with MAP2K6 was reported [9], though not replicated [13]. Other associations for p38 MAPK candidate genes were reported in a sample of 1,102 patients using a generously non-stringent significance threshold of p <0.1 [7].

A major strength of the present study is the use of repeated measures of anti-TNF response at 3 months and 6 months after treatment was started. Previous studies included clinical response from a single time point in standard linear or logistic regression models [59, 11, 12, 1619]. However, assessment of response at a single point in time may not adequately reflect a patient’s response to therapy, as response may fluctuate over time. Hence, using response data from at least two time points is more reliable and clinically relevant. The longitudinal approach enables the use of repeated measures of response from different time-points for each patient, thus increasing the power to detect an association as we have demonstrated, while taking into account within-patient correlation. Patients with missing data at one time point were still included in the analyses, as the GEE uses all available data. However, while the association with the MAP3K7 locus achieved genome-wide significance for SNP rs284511 when using anti-TNF response at only 6 months, the lack of association with this SNP at 3 months led to the GEE model only detecting a borderline significant association.

We did not identify any other overlap between our results and previous findings, possibly due to differences in ethnicity, response variable, i.e., two time points vs a single time point, or duration of anti-TNF treatment (3–12 months in previous studies). Further, the lack of consistent findings between previous studies may also have been the result of differences in a number of factors including study design, clinical outcomes examined (DAS28ESR vs DAS28CRP), specific anti-TNF medications used, concomitant disease modifying anti-rheumatic drugs (DMARDs), sex ratios, and small sample sizes in some cases [9, 12, 19]. Another important difference between studies is the heterogeneity in phenotypes introduced by differences in time from baseline to assessment of response ranging from 3 to 12 months in different studies. As seen in our data, treatment duration is significantly associated with response, and should be given due consideration when using multiple datasets for combined analyses to minimize phenotypic heterogeneity or when comparing results between studies. Further, the low significance thresholds used to identify previously reported associations may have led to false positive associations being included. We also cannot exclude the possibility that our findings may include false positives until they can be replicated in independent datasets or that our sample size was not adequately powered to detect some of the previously reported findings.

The present study has a number of limitations. First, the sample size of 444 patients is modest compared to previous GWAS, which combined data from different populations of European ancestry to achieve large sample sizes [16, 18]. Nonetheless, it represents the largest reported sample size for studies of anti-TNF response among Japanese or any East Asian RA population [18, 24]. Given the lack of pharmacogenomics studies examining associations with anti-TNF response among Japanese RA patients, these results represent an important contribution to the field. Although power may be limited due to the modest sample size, this may in part be compensated for by the ethnic homogeneity of the patient population and the use of response data from two time points as described earlier. Second, the response phenotype in our study included response to three different anti-TNF agents, which may have introduced some bias in responder status because a patient who was a non-responder to one drug, might have responded well to another drug. For example, of the 18 patients who had switched to a new anti-TNF drug at baseline due to inefficacy of the previous drug, 11 had a good response to the second line of anti-TNF agent. However, all patients starting a new anti-TNF agent at study baseline - i.e., those who switched and those who were anti-TNF naïve - were followed for 6 months to assess response as is routinely done in clinical practice. Hence, none of them were switched to a new agent for the duration of the study. To mitigate misclassification bias in the response phenotype that might have arisen from inclusion of three different anti-TNF agents in the analysis, we adjusted for the type of anti-TNF agent used. Third, as there is no gold standard measure to evaluate treatment response in RA, we used the DAS28CRP3 as a surrogate to assess disease activity although it may not be a perfect measure of response. The possibility that such a complex phenotype may be associated only with modest genetic effects has been raised [4, 16]. A comparison of the component variables, i.e., tender and swollen joint counts and CRP levels, showed similar trajectories from baseline to 6 months to the composite DAS28 score. This suggests that no major differences in associations would be expected by focusing on components of the DAS28 as the outcome in our dataset. In order to more closely capture variations in patient response over time, we chose to use repeated measures of the change in DAS28 from baseline. We did not categorize the outcome into European League Against Rheumatism (EULAR) responses as this would have led to a reduction in power. Last, the association between a genetic predictor and clinical response could be confounded by factors that influence response. As far as possible, we adjusted for likely confounders associated with change in DAS28 in our dataset. We did not, however, adjust for variations in drug dosage. We had previously reported that response to anti-TNF therapy is influenced by sex in the long term [60], but that these sex differences were not observed during the first 6 months of treatment. Thus, the lack of an association with sex in the present dataset, which was followed for only 6 months, is in agreement with our previous findings.

Conclusions

In summary, we have identified three chromosomal regions demonstrating significant or borderline significant association with response to anti-TNF therapy among Japanese RA patients. The replication in our Japanese dataset of a previous association with the GFRA1 locus among Caucasians provides evidence for a trans-ethnic association of this locus with anti-TNF response. We have also demonstrated the importance of including response at more than one time point in order to enhance power to detect associations in such pharmacogenomics studies of RA.

Notes

Abbreviations

ΔDAS-3 and ΔDAS-6: 

change in Disease Activity Scores based on 28 joint counts (from baseline) at 3 months and 6 months

ACPA: 

anti-citrullinated protein antibodies

ACR: 

American College of Rheumatology

ADA: 

adalimumab

BACH2: 

basic leucine zipper transcription factor 2

bp: 

base pairs

CRP: 

C-reactive protein

DAS28: 

Disease Activity Scores based on 28 joint counts

DMARD: 

disease modifying anti-rheumatic drug

ETN: 

etanercept

GDNF: 

glial cell line-derived neurotrophic factor

GEE: 

generalized estimating equations

GFRA1: 

glial cell line-derived neurotrophic factor family receptor alpha 1

GWAS: 

genome-wide association study

HWE: 

Hardy Weinberg equilibrium

IL-17: 

interleukin-17

INF: 

infliximab

LD: 

linkage disequilibrium

MAF: 

minor allele frequency

MAP3K7: 

mitogen-activated protein kinase kinase kinase 7

PCA: 

principal component analysis

RA: 

rheumatoid arthritis

RAPA: 

rheumatoid arthritis particle agglutination

RF: 

rheumatoid factor

SJC28: 

swollen joint counts for 28 joints

SNP: 

single nucleotide polymorphism

TAK1: 

transforming growth factor beta-activated kinase 1

TJC28: 

tender joint counts for 28 joints

TNF: 

tumor necrosis factor

UTR: 

untranslated region

WDR27: 

WD repeat-containing protein 27

Declarations

Acknowledgements

We thank the patients for their participation, which made this study possible. This work was supported in part by funds from Chugai Pharmaceutical, Mitsubishi-Tanabe Pharma, Eisai, and Takeda Pharmaceutical (Japan). These funders did not have any role in conducting this study or in interpretation and reporting of results. Dr. Honne was supported by the Jichi Medical University Young Investigator Award.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Division of Rheumatology and Clinical Immunology, Jichi Medical University
(2)
23andme
(3)
Research Unit of Gynecology and Obstetrics, Institute of Clinical Research, University of Southern Denmark
(4)
Department of Obstertics & Gynecology, Odense University Hospital
(5)
Department of Radiology, University of Pennsylvania School of Medicine
(6)
School of Public Health, University of California
(7)
Inoue Hospital
(8)
Center for Genetics, Children’s Hospital Oakland Research Institute

References

  1. Elliott MJ, Maini RN, Feldmann M, Kalden JR, Antoni C, Smolen JS, et al. Randomised double-blind comparison of chimeric monoclonal antibody to tumour necrosis factor alpha (cA2) versus placebo in rheumatoid arthritis. Lancet. 1994;344(8930):1105–10.PubMedView ArticleGoogle Scholar
  2. Moreland LW, Baumgartner SW, Schiff MH, Tindall EA, Fleischmann RM, Weaver AL, et al. Treatment of rheumatoid arthritis with a recombinant human tumor necrosis factor receptor (p75)-Fc fusion protein. N Engl J Med. 1997;337(3):141–7.PubMedView ArticleGoogle Scholar
  3. Barrera P, van der Maas A, van Ede AE, Kiemeney BA, Laan RF, van de Putte LB, et al. Drug survival, efficacy and toxicity of monotherapy with a fully human anti-tumour necrosis factor-alpha antibody compared with methotrexate in long-standing rheumatoid arthritis. Rheumatology (Oxford). 2002;41(4):430–9.View ArticleGoogle Scholar
  4. Umicevic Mirkov M, Janss L, Vermeulen SH, van de Laar MA, van Riel PL, Guchelaar HJ, et al. Estimation of heritability of different outcomes for genetic studies of TNFi response in patients with rheumatoid arthritis. Ann Rheum Dis. 2015;74(12):2183–7.PubMedView ArticleGoogle Scholar
  5. Cui J, Saevarsdottir S, Thomson B, Padyukov L, van der Helm-van Mil AH, Nititham J, et al. Rheumatoid arthritis risk allele PTPRC is also associated with response to anti-tumor necrosis factor alpha therapy. Arthritis Rheum. 2010;62(7):1849–61.PubMedPubMed CentralGoogle Scholar
  6. Potter C, Cordell HJ, Barton A, Daly AK, Hyrich KL, Mann DA, et al. Association between anti-tumour necrosis factor treatment response and genetic variants within the TLR and NF{kappa}B signalling pathways. Ann Rheum Dis. 2010;69(7):1315–20.PubMedView ArticleGoogle Scholar
  7. Coulthard LR, Taylor JC, Eyre S, Robinson JI, Wilson AG, Isaacs JD, et al. Genetic variants within the MAP kinase signalling network and anti-TNF treatment response in rheumatoid arthritis patients. Ann Rheum Dis. 2011;70(1):98–103.PubMedView ArticleGoogle Scholar
  8. Pappas DA, Oh C, Plenge RM, Kremer JM, Greenberg JD. Association of rheumatoid arthritis risk alleles with response to anti-TNF biologics: results from the CORRONA registry and meta-analysis. Inflammation. 2013;36(2):279–84.PubMedView ArticleGoogle Scholar
  9. Krintel SB, Palermo G, Johansen JS, Germer S, Essioux L, Benayed R, et al. Investigation of single nucleotide polymorphisms and biological pathways associated with response to TNFalpha inhibitors in patients with rheumatoid arthritis. Pharmacogenet Genomics. 2012;22(8):577–89.PubMedView ArticleGoogle Scholar
  10. Pavy S, Toonen EJ, Miceli-Richard C, Barrera P, van Riel PL, Criswell LA, et al. Tumour necrosis factor alpha -308G- > A polymorphism is not associated with response to TNFalpha blockers in Caucasian patients with rheumatoid arthritis: systematic review and meta-analysis. Ann Rheum Dis. 2010;69(6):1022–8.PubMedView ArticleGoogle Scholar
  11. Plant D, Prajapati R, Hyrich KL, Morgan AW, Wilson AG, Isaacs JD, et al. Replication of association of the PTPRC gene with response to anti-tumor necrosis factor therapy in a large UK cohort. Arthritis Rheum. 2012;64(3):665–70.PubMedPubMed CentralView ArticleGoogle Scholar
  12. Suarez-Gestal M, Perez-Pampin E, Calaza M, Gomez-Reino JJ, Gonzalez A. Lack of replication of genetic predictors for the rheumatoid arthritis response to anti-TNF treatments: a prospective case-only study. Arthritis Res Ther. 2010;12(2):R72.PubMedPubMed CentralView ArticleGoogle Scholar
  13. Acosta-Colman I, Palau N, Tornero J, Fernandez-Nebro A, Blanco F, Gonzalez-Alvaro I, et al. GWAS replication study confirms the association of PDE3A-SLCO1C1 with anti-TNF therapy response in rheumatoid arthritis. Pharmacogenomics. 2013;14(7):727–34.PubMedView ArticleGoogle Scholar
  14. Zervou MI, Myrthianou E, Flouri I, Plant D, Chlouverakis G, Castro-Giner F, et al. Lack of association of variants previously associated with anti-TNF medication response in rheumatoid arthritis patients: results from a homogeneous Greek population. PLoS One. 2013;8(9), e74375.PubMedPubMed CentralView ArticleGoogle Scholar
  15. Potter C, Gibbons LJ, Bowes JD, Cordell HJ, Hyrich K, Isaacs JD, et al. Polymorphisms spanning the TNFR2 and TACE genes do not contribute towards variable anti-TNF treatment response. Pharmacogenet Genomics. 2010;20(5):338–41.PubMedView ArticleGoogle Scholar
  16. Plant D, Bowes J, Potter C, Hyrich KL, Morgan AW, Wilson AG, et al. Genome-wide association study of genetic predictors of anti-tumor necrosis factor treatment efficacy in rheumatoid arthritis identifies associations with polymorphisms at seven loci. Arthritis Rheum. 2011;63(3):645–53.PubMedPubMed CentralView ArticleGoogle Scholar
  17. Umicevic Mirkov M, Cui J, Vermeulen SH, Stahl EA, Toonen EJ, Makkinje RR, et al. Genome-wide association analysis of anti-TNF drug response in patients with rheumatoid arthritis. Ann Rheum Dis. 2013;72(8):1375–81.PubMedView ArticleGoogle Scholar
  18. Cui J, Stahl EA, Saevarsdottir S, Miceli C, Diogo D, Trynka G, et al. Genome-wide association study and gene expression analysis identifies CD84 as a predictor of response to etanercept therapy in rheumatoid arthritis. PLoS Genet. 2013;9(3), e1003394.PubMedPubMed CentralView ArticleGoogle Scholar
  19. Liu C, Batliwalla F, Li W, Lee A, Roubenoff R, Beckman E, et al. Genome-wide association scan identifies candidate polymorphisms associated with differential response to anti-TNF treatment in rheumatoid arthritis. Mol Med. 2008;14(9–10):575–81.PubMedPubMed CentralGoogle Scholar
  20. Canhão H, Rodrigues A, Santos MJ, Carmona-Fernandes D, Bettencourt BF, Cui J, et al. TRAF1/C5 but not PTPRC variants are potential predictors of rheumatoid arthritis response to anti-tumor necrosis factor therapy. BioMed Res Int. 2015;490295. doi:10.1155/2015/490295. Epub 2015 Mar 5.
  21. Swierkot J, Bogunia-Kubik K, Nowak B, Bialowas K, Korman L, Gebura K, et al. Analysis of associations between polymorphisms within genes coding for tumour necrosis factor (TNF)-alpha and TNF receptors and responsiveness to TNF-alpha blockers in patients with rheumatoid arthritis. Joint Bone Spine. 2015;82(2):94–9.PubMedView ArticleGoogle Scholar
  22. Marquez A, Ferreiro-Iglesias A, Davila-Fajardo CL, Montes A, Pascual-Salcedo D, Perez-Pampin E, et al. Lack of validation of genetic variants associated with anti-tumor necrosis factor therapy response in rheumatoid arthritis: a genome-wide association study replication and meta-analysis. Arthritis Res Ther. 2014;16(2):R66.PubMedPubMed CentralView ArticleGoogle Scholar
  23. Julia A, Fernandez-Nebro A, Blanco F, Ortiz A, Canete JD, Maymo J, et al. A genome-wide association study identifies a new locus associated with the response to anti-TNF therapy in rheumatoid arthritis. Pharmacogenomics J. 2015. doi:10.1038/tpj.2015.31. Epub 2015 Apr 21.
  24. Nishimoto T, Seta N, Anan R, Yamamoto T, Kaneko Y, Takeuchi T, et al. A single nucleotide polymorphism of TRAF1 predicts the clinical response to anti-TNF treatment in Japanese patients with rheumatoid arthritis. Clin Exp Rheumatol. 2014;32(2):211–7.PubMedGoogle Scholar
  25. Yamada R, Yamamoto K. Mechanisms of disease: genetics of rheumatoid arthritis–ethnic differences in disease-associated genes. Nat Clin Pract Rheumatol. 2007;3(11):644–50.PubMedView ArticleGoogle Scholar
  26. Kochi Y, Suzuki A, Yamada R, Yamamoto K. Ethnogenetic heterogeneity of rheumatoid arthritis-implications for pathogenesis. Nat Rev Rheumatol. 2010;6(5):290–5.PubMedView ArticleGoogle Scholar
  27. Okada Y, Terao C, Ikari K, Kochi Y, Ohmura K, Suzuki A, et al. Meta-analysis identifies nine new loci associated with rheumatoid arthritis in the Japanese population. Nat Genet. 2012;44(5):511–6.PubMedView ArticleGoogle Scholar
  28. Viatte S, Plant D, Raychaudhuri S. Genetics and epigenetics of rheumatoid arthritis. Nat Rev Rheumatol. 2013;9(3):141–53.PubMedPubMed CentralView ArticleGoogle Scholar
  29. Okada Y, Wu D, Trynka G, Raj T, Terao C, Ikari K, et al. Genetics of rheumatoid arthritis contributes to biology and drug discovery. Nature. 2014;506(7488):376–81.PubMedPubMed CentralView ArticleGoogle Scholar
  30. Koike T, Harigai M, Ishiguro N, Inokuma S, Takei S, Takeuchi T, et al. Safety and effectiveness of adalimumab in Japanese rheumatoid arthritis patients: postmarketing surveillance report of the first 3,000 patients. Mod Rheumatol. 2012;22(4):498–508.PubMedView ArticleGoogle Scholar
  31. Koike T, Harigai M, Inokuma S, Inoue K, Ishiguro N, Ryu J, et al. Postmarketing surveillance of the safety and effectiveness of Etanercept in Japan. J Rheumatol. 2009;36(5):898–906.PubMedView ArticleGoogle Scholar
  32. Avila-Pedretti G, Tornero J, Fernandez-Nebro A, Blanco F, Gonzalez-Alvaro I, Canete JD, et al. Variation at FCGR2A and functionally related genes is associated with the response to anti-TNF therapy in rheumatoid arthritis. PLoS One. 2015;10(4), e0122088.PubMedPubMed CentralView ArticleGoogle Scholar
  33. Ferreiro-Iglesias A, Montes A, Perez-Pampin E, Canete JD, Raya E, Magro-Checa C, et al. Replication of PTPRC as genetic biomarker of response to TNF inhibitors in patients with rheumatoid arthritis. Pharmacogenomics J. 2015. doi:10.1038/tpj.2015.29. Epub 21 Apr.
  34. Arnett FC, Edworthy SM, Bloch DA, McShane DJ, Fries JF, Cooper NS, et al. The American Rheumatism Association 1987 revised criteria for the classification of rheumatoid arthritis. Arthritis Rheum. 1988;31(3):315–24.PubMedView ArticleGoogle Scholar
  35. Prevoo ML, van 't Hof MA, Kuper HH, van Leeuwen MA, van de Putte LB, van Riel PL. Modified disease activity scores that include twenty-eight-joint counts. Development and validation in a prospective longitudinal study of patients with rheumatoid arthritis. Arthritis Rheum. 1995;38(1):44–8.Google Scholar
  36. DAS-score.nl - Home of the DAS. http://www.das-score.nl/.
  37. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a toolset for whole-genome association and population-based linkage analysis. Am J Hum Genet. 2007;81(3):559–75.Google Scholar
  38. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet. 2006;38(8):904–9.PubMedView ArticleGoogle Scholar
  39. Patterson N, Price AL, Reich D. Population structure and eigenanalysis. PLoS Genet. 2006;2(12), e190.PubMedPubMed CentralView ArticleGoogle Scholar
  40. Delaneau O, Zagury JF, Marchini J. Improved whole-chromosome phasing for disease and population genetic studies. Nat Methods. 2013;10(1):5–6.PubMedView ArticleGoogle Scholar
  41. Howie BN, Donnelly P, Marchini J. A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet. 2009;5(6), e1000529.PubMedPubMed CentralView ArticleGoogle Scholar
  42. Howie B, Marchini J, Stephens M. Genotype imputation with thousands of genomes. G3 (Bethesda). 2011;1(6):457–70.View ArticleGoogle Scholar
  43. Diggle PJ, Heagerty P, Liang K-Y, Zeger SL. Analysis of Longitudinal Data. 2nd ed. Oxford: Oxford University Press; 2013.Google Scholar
  44. Takaesu G, Surabhi RM, Park KJ, Ninomiya-Tsuji J, Matsumoto K, Gaynor RB. TAK1 is critical for IkappaB kinase-mediated activation of the NF-kappaB pathway. J Mol Biol. 2003;326(1):105–15.PubMedView ArticleGoogle Scholar
  45. Sakurai H. Targeting of TAK1 in inflammatory disorders and cancer. Trends Pharmacol Sci. 2012;33(10):522–30.PubMedView ArticleGoogle Scholar
  46. Hammaker DR, Boyle DL, Inoue T, Firestein GS. Regulation of the JNK pathway by TGF-beta activated kinase 1 in rheumatoid arthritis synoviocytes. Arthritis Res Ther. 2007;9(3):R57.PubMedPubMed CentralView ArticleGoogle Scholar
  47. Kooloos WM, Wessels JA, van der Straaten T, Huizinga TW, Guchelaar HJ. Criteria for the selection of single nucleotide polymorphisms in pathway pharmacogenetics: TNF inhibitors as a case study. Drug Discov Today. 2009;14(17–18):837–44.PubMedView ArticleGoogle Scholar
  48. Hammaker D, Firestein GS. "Go upstream, young man": lessons learned from the p38 saga. Ann Rheum Dis. 2010;69 Suppl 1:i77–82.PubMedPubMed CentralView ArticleGoogle Scholar
  49. Clark AR, Dean JL. The p38 MAPK pathway in rheumatoid arthritis: a sideways look. Open Rheumatol J. 2012;6:209–19.PubMedPubMed CentralView ArticleGoogle Scholar
  50. Roychoudhuri R, Hirahara K, Mousavi K, Clever D, Klebanoff CA, Bonelli M, et al. BACH2 represses effector programs to stabilize T(reg)-mediated immune homeostasis. Nature. 2013;498(7455):506–10.PubMedView ArticleGoogle Scholar
  51. Vahedi G, Kanno Y, Furumoto Y, Jiang K, Parker SC, Erdos MR, et al. Super-enhancers delineate disease-associated regulatory nodes in T cells. Nature. 2015;520(7548):558–62.PubMedPubMed CentralView ArticleGoogle Scholar
  52. McAllister K, Yarwood A, Bowes J, Orozco G, Viatte S, Diogo D, et al. Identification of BACH2 and RAD51B as rheumatoid arthritis susceptibility loci in a meta-analysis of genome-wide data. Arthritis Rheum. 2013;65(12):3058–62.PubMedPubMed CentralView ArticleGoogle Scholar
  53. Franke A, McGovern DP, Barrett JC, Wang K, Radford-Smith GL, Ahmad T, et al. Genome-wide meta-analysis increases to 71 the number of confirmed Crohn's disease susceptibility loci. Nat Genet. 2010;42(12):1118–25.PubMedPubMed CentralView ArticleGoogle Scholar
  54. Christodoulou K, Wiskin AE, Gibson J, Tapper W, Willis C, Afzal NA, et al. Next generation exome sequencing of paediatric inflammatory bowel disease patients identifies rare and novel variants in candidate genes. Gut. 2013;62(7):977–84.PubMedPubMed CentralView ArticleGoogle Scholar
  55. Dubois PC, Trynka G, Franke L, Hunt KA, Romanos J, Curtotti A, et al. Multiple common variants for celiac disease influencing immune gene expression. Nat Genet. 2010;42(4):295–302.PubMedPubMed CentralView ArticleGoogle Scholar
  56. International Multiple Sclerosis Genetics C, Wellcome Trust Case Control C, Sawcer S, Hellenthal G, Pirinen M, Spencer CC, et al. Genetic risk and a primary role for cell-mediated immune mechanisms in multiple sclerosis. Nature. 2011;476(7359):214–9.View ArticleGoogle Scholar
  57. Cooper JD, Smyth DJ, Smiles AM, Plagnol V, Walker NM, Allen JE, et al. Meta-analysis of genome-wide association study data identifies additional type 1 diabetes risk loci. Nat Genet. 2008;40(12):1399–401.PubMedPubMed CentralView ArticleGoogle Scholar
  58. Kuno R, Yoshida Y, Nitta A, Nabeshima T, Wang J, Sonobe Y, et al. The role of TNF-alpha and its receptors in the production of NGF and GDNF by astrocytes. Brain Res. 2006;1116(1):12–8.PubMedView ArticleGoogle Scholar
  59. Bian F, Qi H, Ma P, Zhang L, Yoon KC, Pflugfelder SC, et al. An immunoprotective privilege of corneal epithelial stem cells against Th17 inflammatory stress by producing glial cell-derived neurotrophic factor. Stem Cells. 2010;28(12):2172–81.PubMedPubMed CentralView ArticleGoogle Scholar
  60. Jawaheer D, Olsen J, Hetland ML. Sex differences in response to anti-tumor necrosis factor therapy in early and established rheumatoid arthritis – results from the DANBIO registry. J Rheumatol. 2012;39(1):46–53.PubMedPubMed CentralView ArticleGoogle Scholar

Copyright

© Honne et al. 2016

Advertisement