Study population
All RA patients recruited for this study satisfied the American College of Rheumatology diagnostic criteria for RA [6] and had >2 years of follow-up since diagnosis. Also, all patients were Caucasian with all four grandparents born in Spain. Two cohorts of patients from different hospitals in Spain were collected to identify the polymorphisms associated with joint destruction (i.e., discovery phase cohort) and to subsequently validate these SNPs (i.e., validation phase cohort).
This study was undertaken in compliance with the Helsinki Declaration. Informed consent was obtained from all participants, and protocols were reviewed and approved by local institutional review boards. Ethical approval was obtained from the Vall d´Hebron Hospital Ethics Committee.
Discovery and validation phase cohorts
In the discovery phase, 527 patients were recruited from five hospital centers: Hospital Universitario de Asturias, Hospital Clínic i Provincial de Barcelona, Hospital Universitari Vall d'Hebron, Hospital Universitario de Guadalajara, and Hospital del Mar. In the replication phase, RA patients were collected from six different hospitals (n = 705): Hospital de San Rafael, Hospital Universitario La Princesa, INIBIC-Hospital Universitario A Coruña, Hospital Clínico San Carlos, Hospital Universitari Germans Trias i Pujol, and Hospital Regional Universitario de Málaga.
Joint damage scoring
Clinical and epidemiological variables were collected from all patients. Hand and feet radiographic images were obtained from all patients during the inclusion period. Importantly, joint damage was quantified using the same method in all participating rheumatology departments. This S-score method is a joint damage scoring system that captures the level of joint destruction in each patient using a more practical approach compared with other, more complex methods [7, 8]. Given the large number of samples required to analyze the genetic association between IL6R and joint damage, this score was designed to simplify the collection of this trait while reducing the time for acquisition of data in the clinical setting. In the S-score method, six different key areas of the hands and feet are evaluated and quantified from radiographs—distal radius (0/1 points), distal ulna (0/1 points), wrist (0/1 points), metacarpophalangeal joints (0–5 points), interphalangeal joints (0–5 points), and metatarsophalangeal joints (0–5 points)—for each side and for both hands and feet. In each of these areas, the presence of an erosion is considered whenever there is evident loss in the cortical region or loss in the normal bone shape (i.e., referred to as presence of remodeling or loss of bone volume), and contributes 1 point to the S-score. For any patient, the maximal S-score will therefore be 36 (i.e., all areas affected) and the minimal score 0 (i.e., no erosions in hands and/or feet).
To confirm the usefulness of the simplified score to capture joint damage, we first used a sample of RA patients to compare its correlation with a reference score. A total of 139 RA patients with >2 years of follow-up since the diagnosis of the disease were collected, and the degree of joint damage was measured using the S-score and the Sharp–van der Heijde Score (SHS) as the reference score [9]. Although it requires expert training and a substantial amount of time, the SHS is clearly one of the most used scores to quantify joint damage in RA. A highly significant positive correlation (α <0.001) between both scores was considered a validation of the usefulness of the S-score to quantify the level of joint destruction in RA patients. In this group of patients we found a highly significant correlation between the S-score and the SHS reference score (P = 5.42 × 10−20, r
2 (95 % confidence interval (CI)) = 0.67 (0.57–0.75)), and therefore it was considered an informative measure of joint damage.
IL6R tagging SNP selection
TagSNP identification for IL6R was performed using the Tagger method implemented in Haploview software (v 4.2; Broad Institute of MIT/Harvard University, Cambridge, Massachussetts, USA). For this objective, the genotype data from IL6R transcribed region (154,377,669–154,441,926 base pairs (bp) in chromosome 1) were analyzed, including 5000 bp flanking regions at the 5′ and 3′ ends of the gene. To identify the best tagging SNPs, high-density genotyping data from the Hapmap Caucasian European population samples generated 1000 Genomes Project (1KG) were used [10]. Within the IL6R coding region, a total of four splice donor/acceptor variants, three stop gaining variants, three frameshift variants, and 87 missense variants have been identified in the IL6R sequence. Most of these variants, however, are rare (<1 % frequency) or have only been identified at the single individual level (94.3 % variants). Haplotype block analysis identified three main haplotype blocks (Fig. 1a), with a relatively high linkage disequilibrium (LD) between them (r
2 >0.7). Pairwise tagging SNPs were selected (minor allele frequency >0.05, pairwise r
2 >0.8; Additional file 1). TagSNPs tagging <2 SNPs were not selected. A total of five tagSNPs were finally selected that tagged 101 SNPs identified in the IL6R locus.
SNP genotyping
Genomic DNA was isolated from venous blood samples obtained from RA patients using the Chemagic Magnetic Separation Module I (PerkinElmer, Waltham, Massachusetts, USA). SNP genotyping was performed using the TaqMan® genotyping platform (Life Technologies, Carlsbad, California, USA). The TaqMan® assays used to genotype IL-6R tagSNPs were: C_____30047_10 (rs4845618), C__26292294_10 (rs4453032), C__45267034_10 (rs4845374), C__11258914_10 (rs6698040), and C__27968911_10 (rs4379670). Thermal cycle conditions were as follows: 50 °C for 2 minutes and 95 °C for 10 minutes, followed by 40 cycles of 92 °C for 15 seconds and 60 °C for 1 minute. All PCR and endpoint fluorescent readings were performed using an ABI PRISM 7900 HT sequence detection system (Life Technologies). The genotyping error was estimated by genotyping 20 % of the samples in duplicate (error <1 %).
Statistical analysis
The association between IL6R tagging SNPs and joint damage in RA was tested using linear regression with an additive model. Using the linear modeling framework allowed us to control for relevant covariates such as age, years of disease duration, gender, anti-citrullinated protein antibodies (ACPA) and rheumatoid factor (RF) status. All statistical tests were performed using the R statistical software version 3.0.1 [11]. Multiple test correction was performed using the false discovery rate (FDR) procedure.
The meta-analysis of the association results from the discovery and replication cohorts was performed using METAL software [12]. This method uses an inverse variance strategy to weight the association statistics according to their effect size estimates (i.e., beta coefficients from linear regression) and their standard errors from the discovery and validation cohorts of patients.
Haplotype analysis
IL6R haplotype blocks were defined using the model described by Gabriel et al. [13]. Haplotypes were assigned to each individual using PLINK software (version 1.07; Center for Human Genetic Research (CHGR), Massachusetts General Hospital (MGH), and the Broad Institute of Harvard and MIT, Boston, Massachusetts). Analyses of the haplotypes were performed using the same linear regression analysis approach as for the single-marker association analysis.