Abstract

Systemic Lupus Erythematosus (SLE) is a clinically heterogeneous autoimmune disease with strong genetic and environmental components. Our objective was to replicate 25 recently identified SLE susceptibility genes in two distinct populations (Chinese (CH) and Malays (MA)) from Malaysia. We genotyped 347 SLE cases and 356 controls (CH and MA) using the ImmunoChip array and performed an admixture corrected case-control association analysis. Associated genes were grouped into five immune-related pathways. While CH were largely homogenous, MA had three ancestry components (average 82.3% Asian, 14.5% European, and 3.2% African). Ancestry proportions were significantly different between cases and controls in MA. We identified 22 genes with at least one associated SNP (). The strongest signal was at HLA-DRA (; , ); the strongest non-HLA signal occurred at STAT4 (; , ). Most of these genes were associated with B- and T-cell function and signaling pathways. Our exploratory study using high-density fine-mapping suggests that most of the established SLE genes are also associated in the major ethnicities of Malaysia. However, these novel SNPs showed stronger association in these Asian populations than with the SNPs reported in previous studies.

1. Introduction

Systemic Lupus Erythematosus (SLE) is a heterogeneous autoimmune disease, in terms of both clinical presentation and incidence and severity across ethnically diverse populations. Asians are among those with a greater risk of SLE and have more severe disease presentations such as lupus nephritis [1]. SLE has strong and complex genetic components. While several genomewide association studies (GWAS) have been reported for European SLE populations, few Asian GWAS have been performed [24]. Among European identified SLE loci were HLA loci HLA-DRA [5] and ATG5 [5], immune signal transduction loci BANK1 [6], BLK [5], LYN [5], TLR, and IFN pathway related loci IFIH1 [7], STAT4 [8], TNFAIP3 [9], IRF7 [5], IRF8 [7], as well as NCF2 [7], IL10 [10], PHRF1 [5], CD44 [11], ICAM1_ICAM4 [7], TYK2 [7], and UBE2L3 [5]. Loci identified through Asian GWAS include ETS1 [12], SLC15A4 [12], IKZF1 [12], RASGRP3 [12], TNFSF4 [12], and TNIP1 [10].

Malaysia has a population of around 28 million with three major ethnic groups (Malays (60.3%), Chinese (22.9%), and Indians (7.1%)). SLE patients and controls from Malaysia offer a unique opportunity to explore the effect of different ancestral backgrounds [13] on SLE genetic architecture. We explored association of SLE-associated loci identified through GWAS in two majority populations, Chinese and Malays. Given that these cohorts may be admixed, we expect that ancestry proportion may influence SLE association. Although previous studies [1318] reported genetic associations with some candidate genes in Malaysians. To our knowledge, this is the first study which assessed SLE susceptibility genes using large scale targeted fine-mapping on Malaysian populations.

Our objective was to replicate and fine-map genetic association in 25 previously reported SLE susceptibility loci and to assess population structure and individuals admixture in the two ethnically distinct Malaysian cohorts.

2. Materials and Methods

2.1. Subjects and Genotyping

We genotyped 347 cases and 356 controls from the two major Malaysian ethnic groups (Malays (MA) and Chinese (CH)) using the Illumina custom designed ImmunoChip array [19] as part of a separate ongoing genetic association project. The ImmunoChip is a dense fine-mapping genotype array that contains ~196,000 SNPs from 184 genes associated with at least one of 12 autoimmune diseases, including SLE. Genotyping was conducted through the Genotyping Core Facility of the Oklahoma Medical Research Foundation (OMRF), Oklahoma City, USA. Subjects were recruited in compliance with the Internal Review Boards of OMRF and the University of Malaya Medical Centre. All SLE cases fulfilled the ACR criteria for SLE classification [20, 21]. Controls were matched by ethnicity and gender. Our CH cohort included 288 cases and 292 controls (187 males and 393 females); MA included 59 cases and 64 controls (48 males and 75 females) (Supplementary Table 1 in Supplementary Material available online at http://dx.doi.org/10.1155/2014/305436).

2.2. Quality Control

Individuals were removed from analysis if they were genetically related to other study subjects (), as estimated through the relatedness coefficient implemented in GCTA, or if they were outliers (mean ± 2 standard deviations) determined by principal component analysis. SNPs were excluded according to the following criteria: poor genotyping clustering, missing genotype rate greater than 90%, Hardy-Weinberg disequilibrium in controls, or minor allele frequency below 0.5% (Supplementary Figure 1). SNP positions were aligned with HG19. The analysis set contained 6,991 SNPs from 25 previously reported genes genotyped on 580 CH and 123 MA unrelated individuals.

2.3. Population Structure

In order to estimate population structure of our cohorts, we selected 14,134 SNPs with very low intermarker linkage disequilibrium (LD, ). This SNP set was enriched by variants with pairwise allele frequency difference >20%. We merged our cohorts with individuals from the 1000 Genomes Project (103 CEU, 100 JPT, and 101 YRI). We estimated the first ten principal components using GCTA [22], as well as the mean and standard deviation for the first three principal components within each cohort (Figure 1). The same dataset was used to estimate individual admixture proportions in ADMIXTURE [23]. We estimated models of admixture using 1 to 7 ancestry components and determined the optimal admixture model by minimizing the cross-validation error using the Bayesian information criterion and the Akaike information criterion. Mean ancestry between cases and controls was compared with a two-tailed t-test.

2.4. Association Analysis

We performed individual SNP case-control association analysis using a chi-square statistic in PLINK [24]. Given the sample size of our cohorts and that this is a replication study, association was considered significant if (alpha = 0.05). We guarded against type 1 error by performing permutation tests (100,000 permutations). Possible influence of admixture was corrected using a logistic regression model in PLINK [24] with the Asian ancestry proportion as a covariate. We used meta-analysis (Fisher’s combined value, four degrees of freedom) to combine association values from both cohorts. For SNPs which were not significant in either cohort or when odds ratios were not in the same direction, no was calculated. All associated SNPs passed the permutation test (results not shown).

The best SNP was selected for each region starting with the most significant combined . We performed epistasis analysis using PLINK [24] and GAIA [25] in order to identify possible gene-gene interactions. We performed a conditional analysis using a logistic regression model (PLINK) for all significant SNPs from STAT4 and HLA-DRA regions. We used the strongest associated SNP from each loci as the initial conditioned SNP to identify additional independent variants.

In order to check for additional sources of stratification, we used mixed models as implemented on EMMAX [26] (Supplementary Table 2).

RegulomeDB [27] and HaploReg [28] were used to identify functional elements overlapping with the selected SNPs.

2.5. Pathway Analysis

We chose to study five main pathways which contained the majority of the 25 target genes. All of these pathways are reported to be involved in SLE pathogenesis [29]. In order to determine if there were overrepresented pathways in these two cohorts we performed a gene set enrichment analysis (GSEA) weighted by the strength of the meta-analysis association using i-GSEA4GWAS [30]. The objective of this mode of GSEA was to identify the possible biological mechanisms that involve associated loci, and to identify candidate causal SNPs that affect the normal function in these pathways. Since we used a small set of loci we looked for pathways that contained at least two reported genes.

2.6. Power Analysis

We estimated the required sample size for detection of additional association signals in our cohorts using the method developed by Hoggart et al. [31] for alpha = 0.05. This method takes into account the effect of admixture on the probability of identifying an associated variant in a mixed population. The parameters included populations with similar characteristics as CH and MA (admixture proportions of 10% for CH and 20% for MA) and a power of detection of 80%.

3. Results

3.1. Population Structure

Based on the proportions of the 1000 Genomes Project populations we estimated optimal population structure for Asian, African, and European ancestries. CH were very homogenous compared to the MA. As expected, mean Asian ancestral proportion was the highest in both populations (; ), followed by European ancestry (; ) and then African ancestry (; ) (Figure 2). There was a significant mean ancestry difference between cases and controls in MA (case/control: 84.5/80.3: ; 3/3.3: 12.5/16.4: ; ) but not in CH.

3.2. Association Analysis

We identified associated SNPs in 20 previously reported genes in either cohort. However, not all associated SNPs were significant in both cohorts. In CH published non-HLA loci SNPs showed significant association with SLE (Supplementary Table 3), including ETS1 (rs1128334, ), IRF8 (rs2280381, ), TNFAIP3 (rs5029939, ), STAT4 (rs3821236, ), and RASGRP3 (rs13385731, ). In MA, IKZF1 (rs4917014, ), RASGRP3 (rs13385731, ), KIAA1542 (rs4963128, ), TNIP1 (rs10036748, ), and IL21R (rs3093301, ) were significantly associated with SLE. For the HLA locus, we replicated association for rs9271366 (HLA-DRB1_HLA-DQA1  , ; ), consistent with the previous report on Malays and Chinese [13].

For all genes SNPs with the strongest association in this study differed from those previously reported. We identified 22 previously reported genes with at least one associated variant; the most significantly associated SNP for each gene was based on Fisher’s combined value. The strongest association was in the HLA region in the vicinity of HLA-DRA (rs6911777, , , and ). The strongest non-HLA association was observed at STAT4 (rs7568275, , , and ). Also Asian identified TNFSF4 (rs10798269, , , and ) and SLC15A4 (rs6486738, , , and ) were replicated. We identified variants in LD with published variants RASGRP3 (rs13425999, , , , and ), TNIP1 (rs3792782, , , , and ), C7orf72-IKZF1 (rs11185603, , , , and ). Even though these variants have a stronger association signal, they can be explained by their published counterparts.

We also identified a variant in ETS1 (rs76404385, , , and ) that was completely independent of published variant (rs1128334 ). European GWAS identified loci IL10 (rs2232360, , and ), BANK1 (rs17031870, , , and ), PRDM1-ATG5 (rs9398065, , , and ), BLK-FAM167A (rs11782375, , , and ), LYN (rs7828258, , , and ), PDHX-CD44 (rs12362140, , , and ), ITGAM (rs12444713, , and ), NCF2 (rs13306575, , , and ), IFIH1 (rs13023380, , and ), TNFAIP3 (rs5029928, , , and ), PHRF1 (rs4963128, , and ), IL21R (rs8060368, , ), IRF8 (rs34912238, , and ), and ICAM1-ICAM4-TYK2 region (rs12975591, , and ) also had a strong combined association with SLE (Table 1).

Notably, the scales of the odds ratio for rs7568275 (STAT4: , ), rs9398065 (PRDM1-ATG5: , ), rs5029928 (TNFAIP3: , ), and rs76404385 (ETS1: , ) were very close to HLA-DRA levels of OR (rs6911777 , ).

Among the aforementioned 22 SNPs, we identified that rs11782375 (FAM167A_BLK) overlaps with an eQTL that potentially affects gene expression [32]. Additionally, rs13425999 (RASGRP3), rs5029928 (TNFAIP3), and rs11185603 (IKZF1) were identified as likely to affect binding by RegulomeDB [27]. These three SNPs contained enhancer and promoter histone marks in multiple cell types (in particular lymphoblastoid cell type GM12787) and also colocated with DNAse binding sites.

We used conditional analysis to identify multiple independent SNPs for each gene. In particular, STAT4 had rs6740131 (, after conditioning) as an additional independent SNP in CH. In the case of HLA, there were two additional independent SNPs in CH (rs2239806, , and after conditioning and rs532098, , and after conditioning).

3.3. Pathway Related Loci

We identified five important pathways involved in SLE pathogenesis which contained at least one of the 25 genes examined in our study. Both B- and T-cell function and signaling pathways had the greatest number of associated variants (Table 2). Neutrophil/monocyte function and signaling had four significantly associated SNPs, whereas TLR and type I IFN signaling pathways each included five genes with significantly associated SNPs. NFκB signaling also contained SNPs significantly associated with SLE. We did not observe any significantly associated SNPs in DNA degradation apoptosis and clearance of cellular debris pathways.

The only significantly enriched pathway was hsa04514 [cell adhesion molecules (CAMs)]. We derived four causal SNPs that potentially explain enrichment of this pathway, where rs2071554 (nonsynonymous, coding (deleterious) in HLA-DOB) and rs1129740 (nonsynonymous, coding in HLA-DQA1) are candidate causal SNPs through their RECEPTOR_ACTIVITY/TRANSMEMBRANE_RECEPTOR_ACT-IVITY; rs8084 (essential splice site and intronic in HLA-DQB) and rs7192 (nonsynonymous, coding in HLA-DRA) were candidate causal SNPs through TRANSMEMBRA-NE_RECEPTOR_ACTIVITY.

3.4. Gene-Gene Interactions

We did not identify any gene-gene interaction between significant SNPs in either cohort.

3.5. Admixture Correction

We determined potential effects of admixture on associated variants within these pathways by adjusting case-control association analysis with admixture proportions. After admixture correction, only two MA SNPs were no longer significantly associated (). All CH SNPs passed the association threshold ().

4. Discussion

In this fine-mapping study we examined two understudied Malaysian populations to replicate previously known SLE genetic associations and to localize the most associated SNP within known SLE genes. Since SLE heterogeneity may be amplified in admixed populations, we adjusted association for admixture (Asian and European). We also categorized associated variants by pathways involvement and identified particular pathways with accumulation of reported associated variants in our Malaysian populations.

We found no effect of admixture on CH, which was not surprising since they are considered a homogenous population. In fact, the ancestry proportions of European and African were very small, and the minor allele frequencies for the top 25 genes were remarkably similar (Figure 3). Correlation between allele frequencies of CH versus CHB () was higher than MA versus CHB () further supporting our conclusion of the similarity between CH and CHB.

We replicated SLE association in RASGP3 [12], STAT4 [8], TNIP1 [10], IKZF1 [7], IL21R [33], ETS1 [12], and IRF8 [7]. It is not surprising that we did not identify more previously reported loci since the majority of loci were identified from studies of European and European American populations. Given the differences in LD structure between European and Asian populations, we identified new SNPs associated with SLE which could be either causal or in LD with the true causal SNPs within gene.

Associated variants were framed within their possible functional roles in immune-related pathways. The most important SLE-associated pathways in these populations were related to the B- and T-cell function and signaling pathways. We also introduced a causality model for gene set enrichment based on HLA SNPs in the cell adhesion molecules pathway (hsa04514). We identified four SNPs with potential functional effects through an eQTL and histone marks.

Although these results are encouraging, our study is limited due to the small sample size of our cohorts. Given the admixture proportions of these cohorts, we have estimated that we would require at least 1000 cases and controls to identify novel genomewide significant variants () with moderate effects (OR > 1.5) for MA and almost double that for CH. Future large scale admixture mapping with the MA will be especially useful to identify novel SLE susceptibility genes. On the other hand, the CH population can be useful for straightforward association mapping for identifying novel genes or localizing the most likely causal variants.

In conclusion, our high-density fine-mapping on SLE targeted genes is one of the first such undertakings in Malaysian populations. Based on our rigorous analysis, we were able to replicate European and Asian SLE-associated loci in both Malaysian Malays and Malaysian Chinese and were able to identify additional variants that might serve as better tag SNPs for causal variants within these cohorts.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

We are grateful to the affected and unaffected individuals who participated in this study. This work was supported by grants from the US National Institutes of Health (R01AR060366, R21AI094377, and R21AI103399).

Supplementary Materials

Supplementary material contains three tables describing sample demographics; replicated SLE SNPs; and mixed model corrected association results for all our loci. Supplementary Figure 1 describes our approach to quality control through a flowchart.

  1. Supplementary Material