Evaluation of 19 Autoimmune Disease-associated Loci with Rheumatoid Arthritis in a Colombian Population: Evidence for Replication and Gene-Gene Interaction

MATERIALS AND METHODS

We examined association between 19 single-nucleotide polymorphisms (SNP) from 16 gene/genomic regions. These loci were chosen from SNP associated with either RA or at least 1 of the 5 other common autoimmune diseases. Of these, 13 were previously reported to be associated with RA (Table 1; MMEL1, CD244, KIAA1109, IL1/IL2, CDK5, C8orf13-BLK, PHF19-TRAF1, TRAF1-C5, DKFZ, KIF5A, SH2B3, C12orf30, CD40), 4 with SLE (C8orf13-BLK, PHF19-TRAF1, TRAF1-C5, ITGAM), 7 with type 1 diabetes (KIAA1109, ADAD1, IFIH1, IL1/IL2, CLEC16A, CD226, SH2B3), 3 with multiple sclerosis (IL1/IL2, CLEC16A, CD226), and 3 with celiac disease (KIAA1109, IL1/IL2, SH2B3).

We assessed association using a cohort of 353 Colombian patients with RA and 368 healthy controls (543 women, 178 men) enrolled from the Corporación para Investigaciones Biológicas in Medellin. All the study subjects were enrolled at the Center for Autoimmune Diseases Research, in Medellin, and were of Spanish ancestry, originating from the northwestern population of Colombia known as the Paisa community. Historical evidence has documented that individuals from the Paisa community are endogenous, are homogenous, and have very little population stratification^9,10. Written informed consent was obtained from all individuals before enrollment. The ethics committees of Corporación para Investigaciones Biológicas and Oklahoma Medical Research Foundation approved our study. All Colombian patients with RA were diagnosed and followed up by a single rheumatologist, according to the revised American College of Rheumatology criteria¹¹. The average age of onset for patients with RA was 39 (± 12) years.

Genotyping was performed using Taqman SNP genotyping assay in both ABI 7300 and 7500 instruments (Applied Biosystems, Foster City, CA, USA). Polymerase chain reactions were carried out in a total reaction volume of 20 μl, containing 15 ng genomic DNA as a template, 10 μl of Taqman genotyping master media, 0.5 μl of 40× assay mix, and ddH₂O up to 20 μl of final volume.

Statistical analyses were performed with PLINK¹² and Haploview¹³. All SNP were tested for Hardy-Weinberg equilibrium (HWE) in cases and controls separately. Significant differences of allele and genotypic distributions were analyzed using 2-tailed chi-squared tests. To assess the robustness of this significance we performed 100,000 permutation tests. P values, OR, and 95% CI were calculated using Fisher’s exact tests and logistic regression analysis. Power calculations were performed with the genetic power calculator implemented in CATS¹⁴ using a multiplicative model under varying disease allele frequencies (5% to 40%) and genotype relative risks (1.1 to 1.8). Population-attributable risk (PAR) was calculated as PAR = (X – 1)/X. Assuming a multiplicative model, X = (1 – f)² + 2f (1 – f) γ + f² γ², where γ is the estimated OR and f is the frequency of risk allele¹⁵.

Gene-gene interactions were assessed using SNPassoc¹⁶. To adjust for multiple testing we used both Bonferroni corrections and false discovery rate (FDR). For 153 pairwise interactions from 18 SNP with a 5% type-1 error, Bonferroni corrections suggested a threshold of p > 0.0003. We used FDR at a level of 5% (q value, a measure of significance) to adjust p values according to the number of hypotheses tested¹⁷. A q value < 0.05 is considered significant.

RESULTS

The minor allele frequency (MAF) for each SNP is shown in Table 1. All SNP were in Hardy-Weinberg equilibrium in cases and controls except the loci at 10p15, rs12251307 (p_HWE = 2.62 × 10⁻⁶ in cases and 1.71 × 10⁻¹⁰ in controls), which was excluded from further analysis. Our power analysis demonstrated that these samples have more than 80% power to detect an OR of 1.4 for a MAF threshold at 20% or more (using a multiplicative model, at the 5% significance level). We performed a model-based association test for each of the 18 SNP.

We found that both rs13277113 [p = 0.0009; OR 1.46 (95% CI 1.16–1.82)] and rs2736340 [p = 0.0001; OR 1.63 (95% CI 1.27–2.10)], located near the 5’ ends of BLK and C8ORF13, were significantly associated. The MAF for rs2736340 (“T”) and rs13277113 (“A”) was 47% and 46% in cases, and 35% and 37% in controls, respectively. Notably, the MAF for rs13277113 and rs2736340 is much higher in Colombian controls than reported in European-American populations in HapMap data (both ∼25% for both alleles). These 2 SNP were in high linkage disequilibrium (LD; r² = 87%) in both cases and controls. Haplotype-based association testing for these SNP showed significant association (global p = 0.00018), with a higher frequency of “TA” haplotype in cases as compared to controls [44% vs 33%; p = 0.0001, OR 1.40 (95% CI 1.12–1.75)]. Using conditional analysis, neither SNP could independently explain the global p value, which was not surprising given the high LD between these SNP.

We also assessed the association at each SNP between patients who were anti-cyclic citrullinated peptide (CCP) antibody-positive (data were available for 126 patients) and controls. Despite the low power (data not shown) for these analyses, rs2736340 was the only SNP that showed significant association (p = 0.04) in this subgroup analysis with anti-CCP antibody-positive patients.

To explore gene-gene interaction, we assessed all pairwise gene-gene interactions using SNPassoc. The 2 notable pairwise interactions are shown under different models in Table 2. First, an interaction was detected between MMEL1 and C8orf13-BLK under codominant (p = 0.006), dominant (p = 0.01), and log-additive (p = 0.0002) models. The second interaction was detected between CD244 and CDK6 under a codominant model (p = 0.007), dominant model (p = 0.09), and log-additive model (p = 0.01). To correct for multiple testing, we evaluated all pairwise interactions by both Bonferroni corrections and FDR. The interaction p value for MMEL1 × C8orf13-BLK under the log-additive model (p = 0.0002) showed borderline significance after correction for multiple testing using both corrections (p value cutoff = 0.00026) and FDR calculation (q value = 0.06). However, the interaction p value for CD244 and CDK6 was not significant after correction for multiple tests.

DISCUSSION

Recent genome-wide association analyses (GWAS) have implicated several genes associated with RA, and some of these genes have also been shown to be associated with other autoimmune diseases, or vice versa. While the identification of the precise pathways involved in susceptibility to autoimmune diseases will clearly require additional time and effort, integration of data from multiple diseases represents the logical next step in discovering similarities and differences among them. Since these regions of shared autoimmunity are consistently reported from European-derived populations, it is important to confirm these known associations in non-European populations under different geocultural environments.

The 2 associated SNP were located in the C8orf13-BLK region. These SNP were first shown to be associated with SLE in a GWAS, and then replicated in Chinese and Japanese⁷ populations. Recently, 2 studies showed association of SNP from C8orf13-BLK region with RA in a European population. First, a study¹⁸ reported an association between rs2736340 and RA in a European population (p = 5.69 × 10⁻⁹, OR 1.19) with a MAF (allele “T”) of 28% in cases vs 24% in controls, and a PAR of 9%. Second, our recent analysis of the data of the North American Rheumatoid Arthritis Consortium (Genetic Analysis Workshop 16 RA) showed significant association of both rs2736340 (p = 1.45 × 10⁻⁵, OR 1.36) and rs13277113 (p = 3.46 × 10⁻⁶, OR 1.39) with RA. The PAR was estimated at about 14% and 15%, respectively². Of note, Ito, et al¹⁹ recently reported a PAR of 9.2% in a white population; for a Japanese population it was 22% for BLK. Another study²⁰ attempted to replicate recently identified SLE-associated variants in RA samples with European origin. Although they did not replicate association with any of the SNP with RA, they showed association between rs13277113 and a subgroup of patients with RA who had Sjögren’s syndrome (p = 5 × 10⁻⁴). Although functional significance for these SNP remains to be established, reduced expression of BLK has been observed with the risk allele, and detailed mechanisms of this reduced expression are not yet known. It is also possible that these 2 SNP in the 5’ end of the BLK gene may bind transcription factors differentially compared to opposite allele and modulate expression.

We also replicated the association of rs763361 from CD226 (18q22.3) in this Colombian population under the dominant model (p = 0.03). The rs763361 SNP has been shown to be associated with multiple autoimmune disorders²¹ (Table 1) and is an interesting candidate for general autoimmunity. None of the remaining 15 loci was replicated in this Colombian sample. There could be several reasons for this. First, lack of replication of other loci could be due to clinical and genetic heterogeneity of patients with RA, or it could be because of a population-specific effect. Second, it is possible that some of these markers are not the actual predisposing variants, but variants that are in high LD with predisposing variants. A varying LD pattern between Colombian and white populations would be expected. Therefore, without the proper haplotype “tagged” SNP, detecting disease-predisposing SNP may not be possible. A detailed, dense genotyping and association analysis is required to resolve these issues. Third, it is also possible that our study lacked sufficient power to detect modest associations. While we may have had sufficient power to detect a larger effect size (i.e., OR ≥ 1.6), our sample size was underpowered to detect association when OR were ≤ 1.2 (Table 1). Indeed, most of the susceptibility genes detected through GWAS have more modest effects around that range and require a larger sample. For example, CD40 was identified as a novel RA susceptibility gene only after the meta-analysis involving 7322 cases and 18,207 controls^3,22.

We detected evidence of interaction between MMEL1-C8orf13 and BLK regions that retained significance after both Bonferroni correction and FDR. We also detected suggestive evidence for interaction between CD244 (rs3766379) and CDK6 (rs42041) under dominant, codominant, and log-additive models (Table 2). This was not significant after correction for multiple testing, although possible interactions could not be ruled out. Both these SNP have been described to be associated with RA^3,23,24. It is important to note that with our sample size we had a modest power to assess gene-gene interaction. For example, using the genetic power calculator for gene-gene interaction (QUANTO), we assessed that with a disease prevalence of about 1%, and level of significance α = 0.0002, using a log-additive model we had about 7% power to detect interaction between CD244 and CDK6 and about 20% power to detect interaction between MMEL1 and C8orf13-BLK. Extensive search of the microRNA database shows rs3890745 of MMEL1 indicate that this intronic SNP gives very high homology (data not shown) with mature miRNA of various mammals and with humans (hsa-miR-544). We could speculate that this SNP-carrying sequence could code a miRNA that differentially interacts in the promoter region of BLK (rs13277113 carrying alternative allelic sequences) and modulates the expression of the BLK gene.

However, we acknowledge that the gene-gene interaction we observed is not a confirmative but provides a suggestive indication only. Therefore, even in this homogenous population, it is necessary to further evaluate these interactions in another population. Nevertheless, these analyses provide the basis for testing interactions in larger sample sizes of other ethnic groups.

We replicated the association with 3 SNP for RA in a Colombian population and reported evidence of possible interaction between 2 susceptibility genes (MMEL1-C8orf13-BLK) for the first time. These genes are now members of the growing list of confirmed genes with relatively common polymorphisms that have been implicated in the predisposition to multiple autoimmune diseases³. Contributions of these genes to disease predisposition in RA clearly indicated common disease mechanisms between RA and other autoimmune disorders. This information may help us to understand the underlying genetic mechanisms of autoimmunity and may help to identify future molecular targets for disease diagnosis and therapy.