Article Text

Download PDFPDF

Extended report
Genetic variants in IL15 associate with progression of joint destruction in rheumatoid arthritis: a multicohort study
  1. R Knevel1,
  2. A Krabben1,
  3. E Brouwer2,
  4. M D Posthumus2,
  5. A G Wilson3,
  6. E Lindqvist4,
  7. T Saxne4,
  8. D de Rooy1,
  9. N Daha1,
  10. M P M van der Linden1,
  11. G Stoeken1,
  12. L van Toorn1,
  13. B Koeleman5,
  14. R Tsonaka6,
  15. A Zhernakoza1,
  16. J J Houwing-Duistermaat6,
  17. R Toes1,
  18. T W J Huizinga1,
  19. A van der Helm-van Mil1
  1. 1Department of Rheumatology, Leiden University Medical Center, Leiden, Netherlands
  2. 2Department of Rheumatology, University Medical Center Groningen, Groningen, Netherlands
  3. 3Department of Musculoskeletal Sciences, University of Sheffield, Sheffield, UK
  4. 4Department of Rheumatology, Lund University, Lund, Sweden
  5. 5Department of Medical Genetics, Complex Genetics Section, Utrecht, Netherlands
  6. 6Department of Medical Statistics, Leiden University Medical Center, Leiden, Netherlands
  1. Correspondence to Rachel Knevel, Leiden University Medical Center, Department of Rheumatology, Leiden 2333ZA, Netherlands; r.knevel{at}


Background Interleukin (IL)-15 levels are increased in serum, synovium and bone marrow of patients with rheumatoid arthritis (RA). IL-15 influences both the innate and the adaptive immune response; its major role is activation and proliferation of T cells. There are also emerging data that IL-15 affects osteoclastogenesis. The authors investigated the association of genetic variants in IL15 with the rate of joint destruction in RA.

Method 1418 patients with 4885 x-ray sets of both hands and feet of four independent data sets were studied. First, explorative analyses were performed on 600 patients with early RA enrolled in the Leiden Early Arthritis Clinic. Twenty-five single-nucleotide polymorphisms (SNPs) tagging IL-15 were tested. Second, SNPs with significant associations in the explorative phase were genotyped in data sets from Groningen, Sheffield and Lund. In each data set, the relative increase of the progression rate per year in the presence of a genotype was assessed. Subsequently, data were summarised in an inverse weighting meta-analysis.

Results Five SNPs were significantly associated with rate of joint destruction in phase 1 and typed in the other data sets. Patients homozygous for rs7667746, rs7665842, rs2322182, rs6821171 and rs4371699 had respectively 0.94-, 1.04-, 1.09-, 1.09- and 1.09-fold rate of joint destruction compared to other patients (p=4.0×10−6, p=3.8×10−4, p=5.0×10−3, p=5.0×10−3 and p=9.4×10−3).

Discussion Independent replication was not obtained, possibly due to insufficient power. Meta-analyses of all data sets combined resulted in significant results for four SNPs (rs7667746, p<0.001; rs7665842, p<0.001; rs4371699, p=0.01; rs6821171, p=0.01). These SNPs were also significant after correction for multiple testing.

Conclusion Genetic variants in IL-15 are associated with progression of joint destruction in RA.

Statistics from


Rheumatoid arthritis (RA) is an autoimmune disorder that affects 0.5–1% of the population and is associated with significant morbidity, disability and costs for society. Radiographic joint destruction reflects the cumulative burden of inflammation and is conceived as an objective measure of RA severity. The degree of joint destruction varies significantly between patients. The processes behind this difference are incompletely understood. Inflammatory markers and autoantibodies are potent risk factors for joint destruction but account for approximately 30% of the variance in joint destruction.1 Recent data indicate that genetic factors influence the severity of joint destruction in RA.2 Hence, to increase the understanding of progression mediating disease processes, it is relevant to study genetic variants that could predispose for a severe outcome of RA.

In this article, we describe a candidate gene study to investigate the association of interleukin (IL)-15 with the rate of joint destruction. IL-15 is expressed primarily by macrophages, as well as by fibroblast-like synoviocytes and endothelial cells. 3 IL-15 produced by synoviocytes has been implicated in the pathogenesis of RA as it serves as T lymphocyte activator and proliferator in the synovial membrane as well as in the bone marrow.4,,6 IL-15 may be implicated in the perpetuation of synovial inflammation by generating a positive feedback in which activated synovial macrophages or fibroblasts induce continuous T cell recruitment.7 Results of several animal studies supported the notion that IL-15 has a role in RA progression. Ferrari-Lacraz et al generated a mutated IL-15 fused to the constant region of a murine IgG2a, which inhibited the IL-15 receptor. This mutated IL-15-IgG prevented the development of collagen-induced arthritis and also blocked disease progression in an established disease model.8 Daily intraperitoneal injection of mutated IL-15-IgG in collagen-induced arthritis showed reduced clinical scores and reduced cartilage erosions relative to controls. Another study demonstrated that in the absence of IL-15 signalling, several converging mechanisms of osteoclastogenesis were inhibited in mice.9 In patients with RA, IL-15 levels are increased in serum, synovial fluid and bone marrow.4 ,6 Additionally, serum levels of IL-15 correlated strongly to disease activity.10

These data led us to hypothesise that genetic variants in IL15 are associated with the severity of joint destruction in RA. We tested this hypothesis using four data sets of European patients with RA with longitudinal radiological data on joint destruction. All data sets included patients who were diagnosed in a period when treatment strategies were less aggressive and less controlled than today. These conservative treatment strategies made these data sets suitable for the present study as the natural course was less inhibited.

Patients and methods

Study population

Four data sets consisting of adult European patients with RA were studied. RA was defined according to the 1987 American College of Rheumatology criteria in all cohorts except for the Lund cohort where the 1958 criteria were used. X-rays of both hands and feet were available in all cohorts (table 1).

Table 1

Characteristics for each cohort

Leiden Early Arthritis Clinic cohort (Leiden-EAC)

This cohort concerned 600 patients with early RA from the western part of the Netherlands who were included in the Leiden-EAC between 1993 and 2006.1 Patients were included at time of diagnosis and yearly followed up. X-rays were taken at baseline and on yearly follow-up visits for 7 years. In total, 2846 x-ray sets of hands and feet were available. All x-rays were chronologically scored by one experienced reader who was unaware of genetic or clinical data using Sharp–van der Heijde scores (SHSs) on hands and feet.11 A total of 499 x-rays (belonging to 60 randomly selected patients with RA) were scored double. The correlation coefficient (ICC) within the reader was 0.91. The treatment of these patients could be divided into three treatment periods. Patients included in 1993–1995 were initially treated with non-steroidal anti-inflammatory drugs, patients included in 1996–1998 were initially treated with chloroquine or salazopyrine and patients included after 1999 were promptly treated with methotrexate or salazopyrine.


The second set of data involved 275 patients with RA from the northern part of the Netherlands who were diagnosed between 1945 and 2001. The follow-up duration after diagnosis was limited to 14 years. The mean number of x-ray sets (hands and feet) per patient was 3.1 (with a maximum of 8 x-rays per patient). The total number of sets of x-rays was 862. The x-rays were scored chronologically by one of two readers using SHS. ICCs within readers were >0.90, and those between readers were 0.96. The development of joint destruction was significantly different for patients included in the 1990s compared to patients included before 1990. This observation is in line with the introduction of treatment with disease-modifying antirheumatic drug (DMARD).


The third set of data concerned 396 patients with RA from the area of Sheffield, UK. Patients with RA with x-rays available were recruited from the rheumatology department of the Royal Hallamshire Hospital in Sheffield between 1999 and 2006 and were diagnosed as having RA between 1938 and 2003.12 Patients with RA were assessed once during their disease course. The mean±SD disease duration at assessment was 15±11 years (range: 3–65 years). X-rays of hands and feet were scored by one reader using a modification to Larsen's score.13 Ten percent of films were scored twice to quantify the intraobserver variation by a weighted κ score, which was 0.83.12


This cohort concerned 183 Swedish patients with early RA who were prospectively followed up yearly for 5 years, of whom 147 had x-rays and DNA available.14 ,15 Patients were recruited from primary care units in the area of Lund during 1985–1989. X-rays of hands and feet were taken at the start of the study and annually for 5 years, resulting in a total of 781 sets of x-rays. X-rays were scored chronologically according to Larsen by one of two readers.16 The ICC between the readers determined on 105 x-rays was 0.94. In the inclusion period, immediate DMARD treatment was not common and only half of the patients used any DMARD at 5 years of follow-up, most commonly chloroquine, D-penicillamine, sodium aurothiomalate and auranofin.14

All patients gave their informed consent, and approval was obtained from the local ethical committee of each study.

SNP selection and genotyping

The region of IL15, located at chromosome 4q31, plus the regions of the upstream- and downstream-situated haplotype blocks were tagged by the algorithm of haploview.17 No coding and amino acid changing single-nucleotide polymorphisms (SNPs) are known for IL15. One SNP, rs12508866, had a significant association with RA susceptibility in the Wellcome Trust Case Control Consortium data set18 and was therefore forced to be included. Pairwise tagging SNPs were selected from the CEPH/CEU hapmap data set (phase 2, release 21, NCBI (National Center for Biotechnology Information) build 35) using haploview software (minor allele frequency>0.05, pairwise r2<0.8). A total of 25 SNPs captured all SNPs on IL15. Multiplex SNP arrays were designed using Illumina Golden Gate platform, according to the protocols recommended by the manufacturer (Illumina, San Diego, California, USA). Two SNPs (rs9884645 and rs4401531) were excluded as they could not be designed in the multiplex SNP array and no good proxy existed. The SNP selection and the linkage disequilibrium (LD) information are depicted in figure 1.

Figure 1

LD structure between 25 tag-SNPs in IL15. The depicted data are from 600 Leiden-EAC Dutch patients with RA. The numbers present the r2 between the SNPs. The colours refer to D′. Two SNPs could not be designed and no good proxy existed (rs9884645, no. 20; rs4401531, no. 24). One SNP failed typing (rs1961720, no. 7). Significant SNPs in the analyses of Leiden-EAC are marked by an arrow. EAC, Early Arthritis Clinic; IL, interleukin; LD, linkage disequilibrium; SNP, single-nucleotide polymorphism.

Software supplied by Illumina was used to automatically identify the genotypes. Each 96-well plate consisted of one positive and one negative control. In all plates, the positive controls indeed tested positive and the negative controls tested negative. Clusters were evaluated and all doubtful calls were checked; after manually evaluating the spectra of each cluster, we accepted, recalled or rejected the genotypes. At least 12% were assessed in duplicate, with an error rate of <1% for all SNPs except rs7667746, which had an error rate of 3.8%. One SNP failed and two SNPs had a success rate of 75%; the remaining SNPs were typed with a success rate of >98% (Supplementary Table 1). None of the SNPs were out of Hardy–Weinberg equilibrium (p<0.001).

SNPs that were significantly associated with joint destruction in the first data set were genotyped in the other three data sets. SNPs were genotyped with multiplex SNP arrays designed with Sequenom iPLEX, according to the protocols recommended by the manufacturer (Sequenom, San Diego, California, USA). Software supplied by the same manufacturer was used to automatically identify the genotypes. Each iPLEX consisted of at least nine positive and nine negative controls, which were indeed tested positive and negative. All doubtful calls were checked manually. DNA samples that still had >30% failed SNPs after manual checking were excluded from analysis (n=31). At least 5% of the genotypes were assessed in duplicate, with an error rate of <1%. The success rates were all >95%. No SNPs were out of Hardy–Weinberg equilibrium.

Statistical analysis

Associations between genotypes and radiographic joint destruction were analysed. Two phases were carried out. First, an explorative analysis was performed in the Leiden-EAC. In this data set, the tagged SNPs were tested in two ways: additively and recessively. In all data sets, the radiological scores were log-transformed to obtain a normal distribution. Since phase 1 was an explorative phase, no correction for multiple testing was applied yet and SNPs with a p value ≤0.05 were studied in phase 2.

For the analyses in the Leiden-EAC, a multivariate normal regression model for longitudinal data was used with radiological score as response variable. This method analyses all repeated measurements at once and takes advantage of the correlation between these measurements. The effect of time was entered as a factor in the model, to properly capture the mean response profile over time. To test for an association with the rate of joint destruction, we conducted an analysis with the SNP and its interaction with time in the model. The effect of time in the interaction term was linear. Since the analyses were performed on the log scale, the resulting coefficient (β) on the original scale indicates how many fold the joint destruction increased per year in the tested genotype compared to reference genotype and increases per year by the power of the years of follow-up. Adjustment variables were entered based on their univariate association with joint destruction. Adjustments were made for age, gender and the described treatment periods.

For the analyses in Groningen and Lund, we used a multivariate normal regression analysis that was similar to the analysis applied in the Leiden-EAC. Adjustment variables were entered based on their univariate association with joint destruction. The Groningen data set was adjusted for age and inclusion before 1990 and after 1990, as proxy for DMARD treatment. The analysis of Lund was adjusted for age only, since gender and treatment were not associated with joint destruction in this data set.

In the Sheffield data set, each patient had x-rays at one time point. The estimated yearly progression rate was calculated to make the scores comparable to the other data sets.19 This was achieved by dividing the total by the number of disease years at time of x. The disease duration at time of x-ray was available for 391 patients. The SNP association was tested in a linear regression analysis with log-transformed estimated yearly progression rate as outcome variable. No adjustments were applied as none of the tested variables was univariately significantly associated with joint destruction. Also, here, the resulting estimate reflects how many fold the rate of joint destruction increases per year in the presence of a minor allele compared to the absence of this allele. Analyses were performed using SPSS version 17.0.

In the present study, the power is the result of the number of patients and the number of measurements per patient studied. As shown previously, the precision of the estimate increases steadily with increasing numbers of x-rays per patient.20 All three data sets studied to verify the results of phase 1 contained (individually and combined) less x-rays than the initial data set. Consequently, the power to replicate findings in each data set individually as well as in the three replication data sets together was limited in comparison to the large amount of x-rays in the discovery data set. Because of differences in study designs, the separate data sets could not be combined in one analysis directly. Therefore, we decided to test the SNPs in each data set separately, taking advantage of the specific data set characteristics, and to subsequently perform a meta-analysis on the results to determine the association of the SNPs with the rate of joint destruction. A fixed-effects meta-analysis22 with inverse variance weighting was performed in Stata, version 10.1. It is debatable whether correction for multiple testing should be applied. However, multiple testing correction using the Benjamini–Hochberg false discovery rate was performed in phase 2 to reduce the chance of having false-positive findings. p Values ≤0.05 after correction for multiple testing were considered significant.23

Haplotype analyses

Haplotypes of IL15 were studied. Haplotype blocks were defined by Gabriel's method.24 Haplotypes were assigned to each individual using PLINK 1.06 requiring a probability >0.8. Analyses of the haplotypes were performed with methods similar to those used for the analyses of the individual SNPs by now testing the presence of a haplotype compared to the absence of the haplotype.


Discovery phase

IL15 was tagged by 25 SNPs, of which 2 could not be designed, 1 failed typing and 2 did not pass quality control. In total, 20SNPs were analysed in the Leiden-EAC. Five SNPs were significantly associated with joint destruction (minor allele): rs2322182(A), rs7667746(G), rs7665842(G), rs4371699(A) and rs6821171(C) (figure 2). The recessive analyses fitted the data best in all five SNPs. Patients homozygous for the minor allele had respectively 1.04-fold (1.01–1.07, p=0.01), 1.09-fold (1.05–1.13, p<0.01), 1.09-fold (1.04–1.15, p<0.01), 1.09-fold (1.02–1.16, p=0.01) and 0.94-fold (0.90–0.98, p=0.01) higher rates of joint destruction per year as compared to the other patients (table 2). A β of 1.09 per year is equal to a 1.83 (=1.097) higher rate of joint destruction over 7 years. Thus, that is, for rs7667746, the estimated rate of joint destruction in patients carrying both minor alleles was 9% higher per year than that of the other patients; over a period of 7 years, this is equal to an 86% (=0.09237) higher rate of joint destruction. In a subanalysis, analyses were also adjusted for anti-CCP, yielding comparable effect sizes (data not shown).

Figure 2

Depicted are the median SHSs during 7 years of follow-up of patients with different genotypes in phase 1 (Leiden-EAC). Over 7 years, patients with twice the minor alleles of rs2322182, rs7667746, rs7665842, rs4371699 and rs6821171 had 1.31 (1.07–1.60, p=0.01), 1.86 (1.43–2.41, p<0.01), 1.88 (1.33–2.65, p<0.01), 1.80 (1.15–2.82, p=0.01) and 0.64 (0.47–0.87, p=0.01) higher rates of joint destruction compared to patients with only one or no minor allele. A β of 1.09 (rs7667746) per year is equal to a 1.86 higher rate of joint destruction over 7 years, which is similar to the 86% increase in rate of joint destruction. EAC, Early Arthritis Clinic; SHS, Sharp–van der Heijde score.

Table 2

Significant associations between IL15 SNPs and rate of joint destruction as obtained in phase 1 in the Leiden-EAC

Phase 2

The five significant SNPs were assessed in the other three data sets. Patient characteristics are presented in table 1. One SNP, rs7665842, could not be designed by Sequenom. Instead, a good proxy, rs6835391 (r2=0.92), was typed. To refer to this proxy, we used rs7665842. Possibly, due to insufficient power, none of the five SNPs were statistically significant when analysed in each of the other three data sets separately or when tested in the three cohorts combined. Nevertheless, for four of the five SNPs, the effect size in the replication data sets went in the same direction as in the initial data set. The five significant SNPs of phase 1 were analysed in all 1418 patients of all cohorts in an inverse variance weighting meta-analysis in a recessive model. Four SNPs were significant in the meta-analyses: rs7667746 (p<0.001), rs7665842 (p<0.001), rs4371699 (p=0.01) and rs6821171 (p=0.01)(figure 3).

Figure 3

Depicted are the results of the analyses of the individual cohorts and the results of the meta-analyses performed on the four SNPs that were significant in phase 2. The effect sizes are the estimated relative progression rates per year for the presence of twice the minor allele compared to patients with only one or no minor allele. The p values in the graphs are uncorrected for multiple testing. The meta-analyses are based on a fixed-effects model, which are applied to genetic studies to test whether there is a statistically significant effect; generalisability of the effect is of less importance. As result of this choice, the effect size of the meta-analyses should be considered carefully. Consequently, this method is less suitable to estimate the effect size overall. Therefore, the estimated effect of the meta-analysis is depicted in gray. Rs7665842 was typed in the EAC but could not be designed in phase 2. Instead, rs6835391 (r2=0.92) was typed in the other three data sets. EAC, Early Arthritis Clinic; SNP, single-nucleotide polymorphism.

To correct for multiple testing, we calculated the Benjamini–Hochberg false discovery rate q values.23 This revealed that all four SNPs were significant; p values after multiple testing correction were as follows: rs7667746 (p<0.01), rs7665842 (p<0.01), rs4371699 (p=0.02) and rs6821171 (p=0.03).

Haplotype analysis

In all data sets, the three SNPs associated with enhanced destruction (rs7667746, rs7665842 and rs4371699) were in close LD (r2=0.46–0.83, D′=0.97–1). To obtain further insight into the associations found, we performed a haplotype analysis. In the Leiden-EAC, two haplotypes with a prevalence >0.1 were found: GGA and AAC (frequencies: 0.19 and 0.67). Analysis of homozygosity for haplotypes against absence resulted in the following for GGA and AAC, respectively: β=1.09 (1.02–1.16, p=0.01) and β=0.97 (0.95–0.998, p=0.032). Haplotype analysis on the whole gene in the Leiden-EAC resulted in two more blocks (figure 1), which were both not associated with joint destruction (data not shown).

In the three additional data sets, the same haplotypes were made with similar frequencies (frequency: GGA 0.18–0.23 and AAC 0.67–0.69). Meta-analysis of the haplotypes of all data sets resulted in β=1.07 (1.02–1.12, p<0.01) for GGA and β=0.99 (0.97–1.01, p=0.16) for AAC. These effect sizes were comparable to the results of the individual SNPs alone. Therefore, haplotype analyses did not result in additional information on the association discovered.


The variance in joint destruction between patients is considerable and is thus far scarcely understood. We performed a candidate gene study to investigate the association of genetic variants with joint destruction. IL15 (4q31) was chosen as candidate gene since there are emerging data that IL-15 plays a role in perpetuation of inflammation and affects osteoclastogenesis. We tested the association of SNPs tagging IL15 with rate of joint destruction in one data set and tested the significant SNPs subsequently in three other data sets. Finally, an analysis combining the radiological data of all 1418 patients was performed. Four SNPs were observed to associate significantly with progression of joint destruction. One SNP, rs6821171, had a protective effect. Three SNPs (rs7667746, rs7665842 and rs4371699) that were physically closely related associated with a deteriorative effect on joint destruction. Analysis of these three SNPs together in a haplotype analysis did not further elucidate the discovered associations of the single SNPs.

The aim was to study the evolution of joint destruction in the most unbiased manner. The present study uniquely combines four data sets of patients who started treatment in a time when treatment was not as aggressive as nowadays. Since some of the data sets covered a period where different treatment regimens were used, analyses were adjusted for treatment if relevant.

Replication data sets are ideally larger than the initial data set, since effect sizes are generally smaller at a replication stage. However, relatively few large prospective data sets exist with both x-rays and DNA available in conventionally treated patients. For each data set, the number of patients and the number x-rays were lower than the cohort used in phase 1. The data sets available for phase 2 were possibly underpowered to individually replicate findings. In the absence of independent replication, there is an increased risk of a false-positive finding. To minimise this risk, we summarised data in a meta-analysis to test whether the association still holds. An inverse variance weighting meta-analysis yielded significant results for all four SNPs. Importantly, the effects of the SNPs went into the same direction in each data set, which supports the validity of the results.

Different methods to score joint destruction were applied: Larsen and SHS. Although these methods differ in methodology, sensitivity and ranges of scores, both methods are linearly correlated with each other.21 ,25 Therefore, as long as the effects of genetics are tested with the same method within each cohort, the cohort with the different scoring can be combined in one weighted meta-analysis. Nonetheless, the absolute scores provided by both methods are difficult to compare. In the present study, the relative increase in progression rate was evaluated; this relative measure has no units and can therefore be compared.

IL-15 plays a role in T cell proliferation and attraction and has a structural homology with IL-2.26 IL-15 has pleiotropic and physiological activities in both the innate and acquired immune responses. IL-15 induces T cell proliferation, activates NK-cells, costimulates immunoglobulin synthesis by B cells and activates monocytes.27 IL-15 is increased in synovial fluid as well as in bone marrow. It is likely to be involved in the perpetuation of inflammation, which consequently may drive progression of joint destruction in RA.

The present data support the notion that genetic variants in IL15 are involved in the severity of joint destruction. Subsequent studies should elucidate whether the presence of these variants results in difference in IL15 expression, in IL15 activity and in other aspects of IL15 biology.

In conclusion, with a candidate gene approach evaluating patients of four different cohorts, we found an association between genetic variants in IL15 and rate of joint destruction in RA.


The work of AHM van der Helm-van Mil was supported by The Netherlands Organization for Health Research and Development. The work of R Knevel was supported by the Dutch Arthritis Association. The research has been funded by The European Community Seventh Framework Program FP7 Health-F2-2008-223404 (Masterswitch).



  • Competing interests None.

  • Ethics approval The local ethical committee of each cohort has approved this study.

  • Provenance and peer review Not commissioned; externally peer reviewed.

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.