Introduction

Systemic lupus erythematosus (SLE) (OMIM 152,700) is a chronic autoimmune disease that affects multiple organs, and disproportionately affects women and individuals of non-European ancestry1. Candidate gene and genome-wide association studies2,3,4 have successfully identified ∼90 SLE risk loci that explain a significant proportion of SLE’s heritability5,6,7,8. These studies have been largely restricted to populations of European ancestry (EA). Yet, much of the heritability of SLE risk remains unexplained in EA populations, and is largely unknown in other ancestries. Here, we report the results of genotyping large samples of individuals of EA, African American (AA) and Hispanic (Amerindian) American ancestry (HA) on the Illumina Infinium Immunochip (196,524 polymorphisms: 718 small insertion deletions, 195,806 single nucleotide polymorphisms (SNPs)), a microarray designed to perform both deep replication and fine mapping of established major autoimmune and inflammatory disease loci9.

This study identifies 58 distinct non-HLA regions in EA, 9 in AA and 16 in HA. Approximately 50% of the associated regions have multiple independent associations. These 58 regions include 24 novel SLE regions reaching genome-wide significance (P<5 × 10−8). Further, these results localize the association signals in established regions and extended associations to additional ancestries (for example, EA to AA or HA). Adjusting for the associated HLA alleles disentangles a complex multigenic effect just outside of the HLA region. The association between SLE and the risk allele genetic load (risk allele count) exhibits an accelerating nonlinear trend, greater than expected if the loci were acting independently on risk. This nonlinear risk relationship leads us to posit a cumulative hit hypothesis for autoimmune disease. Finally, we report both ancestry-dependent and ancestry-independent contributions to SLE risk.

Results

SLE genetic association study

In total, 27,574 SLE cases and controls from three ancestral groups were genotyped and passed quality control for the Immunochip (AA: 2,970 cases, 2,452 controls; EA: 6,748 cases, 11,516 controls; HA: 1,872 cases and 2,016 controls). Altogether, 146,111 SNPs passed quality control analyses in at least one ancestry (AA: 128,385, EA: 120,873, HA: 120,786). Restricting linkage disequilibrium (LD) to r2<0.2 yielded 46,774 uncorrelated SNPs (union across ancestries) for an estimate of the number of independent tests. To minimize ancestry-specific inflation factors, 3, 4 and 2 admixture factors were included as covariates in the logistic regression model for the EA, AA and HA association analyses, respectively (Supplementary Fig. 1). Inflation factors, scaled to 1,000 cases and 1,000 controls, were λAA,1000=1.03, λEA,1000=1.03 and λHA,1000=1.13 (Supplementary Fig. 2). Power analyses are reported in Supplementary Fig. 3.

Single SNP association

Table 1 shows the number of distinct regions (see Methods) within each ancestry that reached three tiers of statistical significance (Tier 1: P<5 × 10−8, Tier 2: 5 × 10−8<P<1 × 10−6 and Tier 3: P>1 × 10−6 and PFDR<0.05) and lists the number of regions with novel SLE associations. The Tier 1 and Tier 2 thresholds are intentionally more stringent than even the conservative Bonferroni method to reduce the Type 1 error rate on this immune-centric genotyping platform. In total, 5, 38 and 7 distinct non-HLA regions met the Tier 1 threshold of significance for the AA, EA and HA cohorts, respectively; and of these Tier 1 associations, 2, 9 and 2 were novel to SLE regardless of ethnicity or to SLE for a specific ethnicity. An additional 4, 20 and 9 distinct non-HLA regions met the Tier 2 threshold (Fig. 1).

Table 1 Number of non-HLA independent regions per significance tier and ancestry (number of novel regions in parentheses*).
Full size table
Figure 1: Genome-wide associations in SLE.
figure 1

Manhattan plots for (a) European ancestry, (b) African American, (c) Hispanic ancestry, and the (d) meta-analysis. Tier 1 associations are labelled with novel regions highlighted in red. Genome-wide significance (5 × 10−8) is indicated on each plot.

Full size image

European ancestry

Statistically, EA had the most power and 58 regions met Tier 1 or Tier 2 thresholds (Supplementary Data 2). Many are novel SLE risk regions, and others are novel for EA (Table 2). More than 50% of these regions had multiple independent SNPs contributing to the association, based on regional stepwise analyses. In total, 223 distinct associations met PFDR<0.01 (Tables 1 and 2, Supplementary Table 2), which included both well-established and novel associations.

Table 2 Novel ancestry-specific non-HLA associated regions.
Full size table

Novel Tier 1 regions of SLE association in EA and the proximal genes include 4p16 (DGKQ), 6p22 (SLC17A4 and LRRC16A), 6q23 (OLIG3-LOC100130476), 8p23 (FAM86B3P), 8q21 (PKIA-ZC2HC1A) and 17q25 (GRB2). Of the 20 EA Tier 2 associated regions, 16 appear novel to SLE.

African American ancestry

The AA sample was powered to detect OR=1.1 to 1.2 at α=1 × 10−6. In addition to known regions in AA, novel AA regions identified include 5q33 (PTTG1-MIR146A), 6p21 (UHRF1BP1-DEF6) and 16q22 (ZFP90) (Tables 1 and 2; Supplementary Data 2). The 8p11 (PLAT) association is novel to SLE and was not observed in HA or EA as it was nearly monomorphic in both populations. The 1q25 region in AA is near the known anti-dsDNA-rs2205960 association between TNFSF4 and LOC100506023 in non-AA samples. The association at rs6681482 (P=8.11 × 10−7, OR=0.73) within LOC100506023 appears independent and separated from the TNFSF4 associations by a recombination hotspot. Three SNPs in this region met the stepwise significance threshold, but the strongest association in EA (rs2205960) was not genome-wide significant in AA (OR=1.35, P=7.39 × 10−4). The association with rs2431697 (OR=0.76, P=1.27 × 10−12) at 5q33 was previously associated with SLE and anti-dsDNA in EA, but not in AA (ref. 10).

Hispanic ancestry

HA samples had comparable power to the AA sample but exhibited more (nine versus four) novel associations at the Tier 1 and Tier 2 thresholds (Tables 1 and 2). Many regions had multiple independent associations, including cases of previously reported regions exhibiting additional novel loci. Novel Tier 1 regions include 14q31 (GALC) and 16p13 (CLEC16A). Novel Tier 2 regions include 3p11 (EPHA3-PROS1), 6p21 (TCP11-SCUBE3), 6q25 (RSPH3), 12q15 (DYRK2-IFNG), 12q21 (SYT1), 16q21 (CSNK2A2-CCDC113) and 22q12 (C1QTNF6). Only the 16p13 locus is associated in AA and EA.

Chromosome X

None of the 442 chromosome X SNPs, predominantly in Xp22 and Xq28, met Tier 1 or Tier 2 thresholds of significance. The strongest evidence of association was in females at Xq28 within GAB3 (Supplementary Fig. 4; rs2664170; EA: OR=0.89, P=0.0009; AA: OR=0.90, P=0.13; HA: OR=0.90, P=0.33; Meta P=1.23 × 10−4).

Two-way interactions among associated SNPs

No SNP–SNP interactions met the Bonferroni threshold (P=1 × 10−9) (see Methods).

Human leukocyte antigen region

SNP analyses within the HLA region provided strong evidence of association with SLE across groups (Fig. 2). These associations are complicated by the region’s extended LD between SNPs and classical HLA alleles. Supplementary Data 3 and Supplementary Fig. 5 summarize the posterior probability distributions for the imputed four-digit HLA alleles in HLA-A, -B, -C, -DQA1, -DQB1, -DPB1 and -DRB1.

Figure 2: HLA SNP associations with and without adjustment of classical HLA alleles.
figure 2

SNPs spanning the extended MHC region showed significant associations across (a) European ancestry, (b) African American, (c) Hispanic ancestry, and the (d) meta-analysis. The classical HLA alleles, from the ethnic-specific stepwise-models (Supplementary Data 5), accounted for a majority of the MHC SNP signals. For each plot, the Tier 1 threshold, P≤5 × 10−8, is indicated by the red line. Associations, downstream in 6p21 spanning UHRF1BP1-DEF6 were largely unaffected after adjusting for classical HLA alleles and appear independent of the MHC.

Full size image

HLA allele associations

HLA allele associations for each ancestry and for multi-ancestral meta-analysis are shown in Supplementary Data 4. To disenable regional LD effects, ancestry-specific stepwise logistic modelling was used to identify the set of alleles with unique HLA contributions to SLE risk (that is, risk or ‘protective’ alleles associated even after adjusting for other SLE-associated HLA alleles) (Supplementary Data 5). To account for HLA alleles contributing even nominal effects, the models’ entry and exit criteria were set to P≤0.01 (see Methods). The final models contained both risk and ‘protective’ alleles. In both the single-allele and multi-locus models, class II alleles exhibited the greatest association with SLE. The DR3 (DRB1*3:01-DQA1*05:01-DQB1*02:01) and DR15 (DRB1*15:01/03-DQA1*01:02-DQB1*06:01) haplotypes had the most significant class II risk alleles across populations.

SNP associations after adjusting for HLA alleles

Only two SNPs showed evidence of association with SLE (Supplementary Data 6) after adjusting for the HLA alleles identified in the stepwise modelling (Fig. 2). Specifically, for EA these SNPs are, rs1150755 (OR=1.33, P=3.10 × 10−8) within TNXB and rs9273448 (OR=0.64, P=2.39 × 10−8) within HLA-DQB1 (Supplementary Data 6 and Supplementary Fig. 6). These associations had comparable ORs in the AA and HA cohorts, except in HA for rs9273448. Transancestral meta-analysis showed stronger association at both loci (Supplementary Data 6 and Supplementary Fig. 6). Whether these residual associations reflect novel loci or imperfect imputation requires additional study.

Compound risk allele heterozygosity

In several autoimmune diseases, including lupus11, having two different risk alleles (compound risk allele heterozygosity) generates greater disease risk than having two copies of the same risk allele12,13. In SLE, there are two primary risk haplotypes (DRB1*3:01-DQA1*05:01-DQB1*02:01 and DRB1*15-DQA1*01:02-DQB1*06:01), which are comprised of alleles in strong linkage disequilibrium. Thus, we selected DRB1*03:01 and DR*15 (DRB1*15:01 in EA & HA; DRB1*15:03 in AA) as tagging alleles to evaluate risk allele heterozygosity. Supplementary Data 7 summarizes the genotypic associations and contrasts the effects of risk allele homozygosity, heterozygosity, and compound heterozygosity. In both EA and AA, compound risk allele heterozygosity (DRB1*03:01/*15 provided greater risk than homozygosity for either individual risk allele (that is, DRB1*03:01/03:01; 15/15); these effects are consistent in direction but not significant in HA. Transancestral meta-analysis strongly supports that the risk for compound heterozygotes is greater than homozygotes for any individual allele (P03:01=1.79 × 10−10; P15:01=4.65 × 10−28). While there was not conclusive evidence of a statistical interaction for people having these two risk alleles in EA (P=0.07), AA (P=0.06), or HA (P=0.50), the lack-of-fit test supported the dominance model of risk (departure from additivity; see Methods) for an individual DR3 (EA P=7.90 × 10−109; AA P=0.06; HA P=5.14 × 10−10) and DR15 (EA P=5.79 × 10−26; AA P=3.99 × 10−13; HA P=3.25 × 10−11) SLE risk alleles.

HLA clustering by amino acid

HLA alleles with high sequence similarity, but contrasting ORs, suggest the potential presence of key amino acids influencing disease risk. As expected, clustering amino acid sequences resulted in most two-digit allele subtypes residing within the same clusters (Fig. 3 and Supplementary Fig. 7). When evaluating SLE associations of the three ancestries across these sequence clusters, several noteworthy patterns emerged.

Figure 3: Clustering of HLA Class II alleles by amino acid sequence similarity.
figure 3

For (a) DRB1, (b) DQA1, and (c) DQB1, the odds ratios for each cohort are superimposed on the cluster if the SLE association P-value was less than 0.01. Alleles that were present in the multi-locus model from the stepwise procedure are also denoted. This process aims to identify clusters with shared SLE risk or not-risk odds ratios across the three cohorts. Such clusters help identify potential amino acid sequences contributing to SLE risk. For example, DRB1*15:01 and 15:03 are clustered amongst protective alleles, suggesting presence of specific amino acids differentiating risk (Supplementary Figs 8 and 9).

Full size image

The two primary DRB1 risk alleles, DR3 and DR15 clustered separately, suggesting comparative amino acid dissimilarity. Notably, the closest-clustered neighbours to each risk allele conferred non-risk in these three ancestries. Multi-sequence alignment distinguished the unique or less common amino acids among risk alleles (Supplementary Figs 8–10). Unique to risk alleles DRB1*15:01 and *15:03 were the amino acids Ser-1 (signal peptide), Phe47 and Ala71. Three-dimensional modelling of DRB1 (Supplementary Fig. 8b,c) reveals that these differences mostly reside within the peptide-binding pocket, creating a space of non-polar (hydrophobic) residues, unlike the polar-residue (hydrophilic) space of Tyr47 and Arg71 or Glu71 provided by non-risk alleles within this cluster (Supplementary Fig. 9). Residue 71, among the most variable residues in DRB1 (ref. 14), has been implicated in other diseases15. Among non-risk alleles with at least 95% identity to DRB1*03:01, the only amino acid unique to this risk allele was Tyr26 (Supplementary Fig. 10). DRB1*03:01 amino acids shared by less than half of the non-risk alleles in this cluster are highlighted in Supplementary Fig. 10 and are concentrated between positions 70–77, spanning the designated ‘Shared Epitope’ region16,17.

One predominant DQA-DQB1 pair of SLE risk alleles exists per evolutionary DQ-sublineage (Fig. 3b,c)18. In the DQ2/3/4 sublineage, DQA1*05:01 confers risk across the three cohorts and its heterodimer counterpart, DQB1*02:01, confers risk in EA and HA, but not significantly in AA. Within the DQ5/6 sublineage, both DQA1*01:02 and DQB1*06:02 yield SLE risk across all three cohorts. Comparison of DQA1*01:02 to its closest-related alleles (Supplementary Fig. 11) reveals that DQA1*01:02 (DR15) uniquely encodes a Met207 versus Val207. DQA1*05:01 encodes a polar Thr13 compared to the non-polar Ala13 found in DQA1*05:05 (DR3) and DQA1*05:03 (Supplementary Fig. 12). Identification of specific risk residues was less distinct for the DQB1 risk alleles.

Gender-HLA and genome-wide SNP-HLA interaction

There was no evidence that the risk of SLE differed by gender at any HLA alleles or of a significant SNP-by-HLA allele interaction anywhere across the genome (PFDR>0.05).

Transancestral mapping and top meta-analysis regions

The three-ancestry meta-analysis identified additional SLE-associated regions and was particularly informative for 22 regions, including 11 novel regions, 3 published regions that now meet genome-significance, a complex multigenic region identified by adjusting for HLA alleles and 7 well-established regions more sharply localized by transancestral mapping or novel to these ancestries (Tables 3 and 4; Supplementary Figs 13–15). Supplementary Data 8 and Supplementary Fig. 16 show additional regions that only met genome-wide significance in the meta-analysis. Supplementary Data 9 lists any region with meta-analysis PFDR<0.001.

Table 3 Novel non-HLA associated regions identified by transancestral meta-analysis.
Full size table
Table 4 Tier 1 non-HLA meta-analysis regions noted for transracial mapping.
Full size table

On 1p31, rs3828069 is within an intron of IL12RB2 (OR=0.85, P=1.77 × 10−9) and has evidence of association in all three ancestries. Although IL12RB2 is implicated in multiple autoimmune diseases19,20, this specific SNP association with SLE is novel. The 2p16 region exhibited a novel SLE association at rs1432296 (OR=1.18, P=1.34 × 10−8) near PAPOLG-LINC01185, which includes REL. A linkage region at 4p16 (ref. 21) contained a strong novel association for rs3733345 (OR=0.89, P=1.83 × 10−11); EA dominated the association, but with significant support from HA and AA. On 8q21, rs4739134 is near PKIA-ZC2HC1A (OR=1.12, P=3.47 × 10−8) and the AA helped localize the association. The region about 16q13 (PLLP-CCL22) exhibited modest association in individual ancestries, but reached genome-wide significance for rs223889 (OR=1.21, P=1.08 × 10−8) in the meta-analysis. Similarly, rs137956 (OR=0.88, P=5.0 × 10−8) on 22q13 between ENTHD1 and GRAP2 was supported across all three ancestries. We bioinformatically explore three additional novel regions.

The meta-analysis about 16q22 (rs1749792; OR=1.14, P=3.66 × 10−11) near ZFP90 had strong support from both EA and AA, with AA samples localizing the association (Supplementary Fig. 13l). While previously identified in a Chinese cohort, this is the first significant association within EA and AA8. Within this region, 27 additional SNPs had a meta-analysis P value within one order of magnitude of the maximum association, rs1749792. These 28 SNPs span an interval of 44.6 kb, narrowed from the 100 kb associated region in EA. RegulomeDB22 and HaploReg4.1 (ref. 23) identified 4 of these SNPs with a RegulomeDB score of 1f and 1 with a RegulomeDB score of 2f, indicating they were eQTLs and transcription factor binding sites. HaploReg4.1 showed these five SNPs were enhancers and promotor histone marks in multiple tissues. Interestingly, one of these five, rs1170445, is in high LD with rs1749792 (R2EA=0.99, R2AA=0.84, R2HA=0.99). Here, the G allele is the risk allele and creates a CpG site in the promoter region. In GTEx, the G allele corresponds to lowest gene expression. Hence, when methylated, this variant should result in decreased gene expression of ZFP90. The rs1170445-ZFP90 expression association was reported in GTEx for whole blood (P=1 × 10−47) and several other tissues (that is, spleen, skeletal muscle, brain cortex, lung, testis and EBV-transformed lymphocytes). Huang et al.24 found expression of ZFP90 in Jurkat T cells led to decreased expression of IL2 and interferon. Furthermore, they found that ZFP90 protein binds to IL2 and interferon gamma promoters.

SLC15A4 was associated with SLE in the EA cohort and localized by the AA signal in the meta-analysis. The top EA signal was supported by a 43.7 kb region of SLE-associated SNPs exhibiting P values within one order of magnitude of the top signal. The meta-analysis narrowed the region of association to four SNPs, spanning 9.5 kb around rs1059312 (Supplementary Fig. 15j). rs1059312 is an eQTL for SLC15A4 and three supporting SNPs (rs2291349, rs4760593 and rs11059916) altered CpG sites. The region has been previously reported in Asian populations25,26; but this is the first instance of genome-wide significance in EA (P<5 × 10−8)26.

On 17q25 near GRB2, rs8072449 (OR=0.84, P=1.19 × 10−11) had modest support in each ancestry, but met genome-wide significance and better localization in the meta-analysis. rs8072449 is an eQTL for GRB2 (Supplementary Fig. 13m). There were eight additional SNPs with a meta-analysis P value within one order of magnitude of the maximum association, and the transancestral analysis reduced the interval of association from 93 to 82 kb. The best RegulomeDB scores for these 9 SNPs was 1f for rs7219, reflecting rs7219 as a known cis-eQTL (NUP85, MIF4GD, MRPS7), a transcription binding site and within a DNase peak; in total 7 of the 9 SNPs were reported in transcription binding sites. Interestingly, the top associated SNP, rs8072449, breaks a CpG site and 6 others either end or begin a CpG site. Hence, 7 of the 9 top associated SNPs make or break a CpG site and several are transcription binding sites. Of the 147,111 Immunochip SNPs that passed quality control analyses, only 30% begin or end a CpG site. Although this is a novel SLE association, GRB2 reportedly regulates SHP2 activity27,28, a potential contributor to SLE pathogenesis29.

A few novel regions, sparsely mapped on the Immunochip, reached genome-wide significance in the meta-analysis and merit further fine-mapping efforts. These include rs6886392 on 5q21 (OR=1.13, P=4.08 × 10−9), rs11788118 on 9q22 (OR=0.88, P=1.53 × 10−8) and rs13344313 on 19p13 (OR=0.90, P=1.07 × 10−8).

Additional loci not previously reported as having genome-wide significance for SLE in these ancestries now do so in the meta-analysis (Table 4). On 4q27, rs11724582 (OR=0.88, P=1.71 × 10−8) is near IL21, a known SLE risk locus30,31. IL21 is up-regulated by oestrogen and is produced by T follicular helper cells which stimulates B-cells to differentiate into autoantibody-secreting cells; however, there was no evidence of a SNP-by-gender interaction in any ancestry (P>0.40). The SNP rs2431098 (OR=1.19, P=3.29 × 10−21) at 5q33 between PTTG1 and MIR146A has an r2=0.52 with rs2431697, a SNP correlated with down-regulation of MIR146A32.

The 6p21 region is potentially confounded with nearby HLA associations. The advantages of using multiple ancestries in this study are exemplified by modelling of SNPs in the 6p21 region where three separate ancestry-specific signals were identified after adjusting for HLA alleles. The results show associations at previously reported UHRF1BP1 and two novel loci within the SCUBE3-DEF6 region (Fig. 2 and Supplementary Fig. 13e,f).

The transancestral meta-analyses of several previously established SLE associations provided important localization, and increased the number of independent signals or novel transancestral effects. These included: 1q25 (TNFSF4-LOC100506023), 1q25 (NMNAT2-SMG7-NCF2), 7q32 (IRF5-TNPO3), 8q12 (LYN-RPS20), 11p13 (PDHX-CD44) and 20q13 (NCOA5-CD40) (Table 4, Supplementary Fig. 15).

Admixture and population frequencies of SLE-associated SNPs

Clustering risk allele frequencies for Tier 1 and 2 SNPs in cases across EA, AA, and HA yielded three groups of SNPs: comparable allele frequencies in all three ancestries (75 SNPS), increased frequency in AA cases (40 SNPs), and reduced frequency in AA cases (66 SNPs) (Fig. 4); the latter two clusters show increased and decreased AA-ancestral contribution, respectively. Higher frequency risk alleles tend to exhibit comparable frequencies across ancestries; the rarest alleles were largely grouped in the reduced AA-ancestral cluster. When comparing admixture averages for risk alleles, AA exhibited the highest deviations from mean admixture estimates and EA, the lowest (Fig. 4; Supplementary Data 10). Deviations from average admixture in risk alleles were significantly weighted to higher proportions of CEU versus YRI in AA (P=8.36 × 10−12) and HA (P=2.44 × 10−4) (Supplementary Data 11), further suggesting increased European ancestry for risk alleles. When aligned to allele frequency information, highest CEU proportion deviations in AA and HA resided in the decreased-AA cluster, while the YRI proportion deviations resided in the increased-AA cluster. Thus, SLE risk alleles with a low frequency in AA are correlated with European admixture. Of the 181 Tier 1 and 2 SNPs, only in two regions were the top associated SNP (rs1804182 AA Tier 1 and rs11845506 HA Tier 2) nearly monomorphic (frequency<0.003) in the other ancestral cohorts. This suggests that most of the ancestry-specific SNP associations were not driven by the presence of monomorphic alleles in the non-discovery cohorts. These allele patterns are further illustrated in Fig. 4.

Figure 4: Ancestral landscape of SLE risk alleles.
figure 4

Clustering by relative allele frequency yields three distinctive categories for SLE risk alleles: comparable frequencies across populations, increased frequencies in AA, and decreased frequencies in AA. The comparable frequency grouping contained the most risk alleles, of which, many were common alleles. This cluster had the smallest deviations from average admixture proportions, across the three cohorts. The increased frequencies in AA alleles exhibited moderate deviations towards greater AA-ancestral contribution. The largest deviations from average admixture were found within alleles exhibiting decreased frequencies in AA. These alleles were enriched for admixture deviations of increased CEU-ancestry. The patterns across relative allele frequencies reveal that ancestry-specific associations are largely not driven by monomorphic SNPs in other populations.

Full size image

Genetic load and SLE risk

To explore effects of the number of risk polymorphisms on SLE risk, we computed the genetic risk allele load (unweighted and β-weighted (β=log(OR)), see Methods). Here, a set of ORs that contrasted the lowest 10% of the risk-allele count distribution with a sliding window of 20 unweighted, or 4 weighted, counts was computed; these logistic models adjusted for admixture. The pattern of the sliding window ORs was different across ancestries (Fig. 5 and Table 5). Specifically, in 2,000 EA cases and 2,000 EA controls that were independent from the discovery set, a strong and nonlinear effect emerged, with ORunweighted>30 and ORweighted>100 for the highest load groups. In fact, there was a nonlinear trend in the log(OR) (that is, β parameter denoting slope) with a greater than additive effect at the highest quarter of the genetic load range (Supplementary Fig. 17); this pattern suggests that the effect of at least a subset of the alleles is greater when the overall genetic load is high. HA and AA showed markedly smaller ORs (between 3 and 10), reflecting the reduced predictive ability of EA-identified SLE risk loci in non-EA populations and the lack of capturing non-EA SLE risk loci on the Immunochip.

Figure 5: The non-additive effect of EA risk-allele genetic load on SLE risk.
figure 5

The cumulative effect of EA SLE-risk alleles (cumulative hits) on an individual's risk of SLE is greater than if the individual SNPs were acting independently/additively. (a) The genetic load was computed as the sum of the number of EA risk variants from the Tier 1, 2 or 3 SNPs that met the region-specific stepwise modelling (see Online Methods). In the AA, HA and an independent set of 2,000 EA cases and 2,000 EA controls, the samples with the lowest 10% in risk-allele counts were identified and formed the baseline comparison group. Using a moving window of 10 in the allele count, the odds ratio for that window relative to the lowest 10% was computed and graphed. (b) The process was repeated for a weighted sum of the number of EA risk-allele variants. Here, the alleles are weighted by the natural logarithm of the odds ratio for that SNP’s association with SLE. The corresponding moving window for the weighted genetic load used a window size of 3. Supplementary Fig. 17 plots the natural logarithm of the odds ratio (instead of the odds ratio) of genetic load versus SLE risk.

Full size image
Table 5 Genetic Load and SLE risk.
Full size table

The total non-HLA weighted genetic load was correlated with an earlier age at SLE diagnosis in EA (rSpearman=−0.14, P=0.0001), and HA (rSpearman=−0.10, P=0.0012), but not AA (rSpearman=0.04, P=0.54). Kaplan–Meier curves in the EA showed separation accelerates at ∼35 years (Supplementary Fig. 18). The HLA-based genetic load was not correlated with age of onset (P>0.05) in any ancestry.

Mapping SNP associations to eQTLs

Many SLE-associated SNPs are, or are in LD with, cis eQTLs (Supplementary Data 12 and Supplementary Figs 13–16) and potentially link associations with specific genes. In ancestry-specific eQTL analyses (Supplementary Data 12), EA yielded 96 unique SNPs or their proxies mapping to 193 unique genes, followed by HA (22 unique SNPs; 34 genes) and AA (10 unique SNPs; 17 genes). eQTL analyses based on the meta-analysis SNPs yielded 107 unique genes, identified by 40 SNPs (or their proxies), mostly from whole blood, monocytes or B-cell derived LCL (Supplementary Data 12). Novel and previously implicated SLE genes were identified (for example, BANK1, IRF5). Interestingly, a number of SNPs were associated with expression levels for multiple genes. For example, four SNPs were associated with expression levels of at least three genes, and one SNP, newly associated in this study (rs8072449; 17q25), were associated with expression levels of eight genes. Thus, some associated SNPs, either directly or via LD with proxy SNPs, contribute to disease by modifying expression levels of multiple genes, potentially through transcription binding sites. Supplementary Data 13 and 14 provide predicted functional characterization of the 206 SNPs from Tiers 1 to 2 that are in RegulomeDB and HaploReg. These predictions are informative for generating hypotheses that can be experimentally tested.

Discussion

Applying the Immunochip to these multi-ancestral SLE case-control samples has identified 24 novel SLE-risk regions, replicated established SLE-risk loci and extended their impact into other ancestries, and refined association signals via transancestral mapping. Over 50% of associated regions had multiple independent SNP associations. Many of these associations were linked via eQTL analysis to specific genes, a process that can accelerate discovery of critical pathways. The contrast of associations and genes across ancestries documents numerous ethnic-specific associations the ancestral diversity in SLE etiology; for example, HA regions not showing equivalent associations in EA include 3p11 (EPHA3-PROS1), 6q25 (RSPH3), 12q15 (DYRK2-IFNG), 12q21 (SYT1), 14q31 (GALC), 16q21 (CSNK2A2-CCDC113) and 22q12 (C1QTNF6). In total, these results underscore the shared and distinct genetic profiles of SLE relative to other autoimmune diseases.

To understand disease biology and prevalence across populations, distinguishing shared versus ancestry-specific associations is important because an allele identified in one population is likely relevant in others33. Clustering by allele frequencies in cases and comparing risk allele admixture estimates, three clusters emerged: (1) alleles with comparable frequencies across populations without strong deviations in average admixture, (2) alleles with increased AA-ancestral contribution and (3) alleles with reduced AA-ancestral contribution and increased CEU admixture. The increased European ancestry observed in less common AA risk alleles likely reflects complex demographic histories and admixture patterns.

The nonlinear nature of how genetic load affects SLE risk leads us to posit the cumulative hit hypothesis for autoimmune diseases. That is, in our current environment the immune system can absorb, with a modest increase in risk, individual risk polymorphisms. But as the number of risk variants increases, the system becomes overwhelmed and immune dysregulation occurs. Currently, it is unclear whether it is the entire genetic load or only a subset of variants driving the nonlinear association. In addition, increasing genetic load correlates with an earlier age of disease onset. These hypotheses are testable within specific and across autoimmune diseases given their shared genetic architecture.

Despite the large sample size, there was no robust evidence for SNP-gender, SNP–SNP or SNP–HLA allele interactions, suggesting that pairwise-interactions among these Immunochip loci are not a major source of missing heritability. While the lack of pairwise interactions across the immune-centric loci may be surprising given the statistical power of the study, the current analysis does not preclude higher-order interactions; albeit agnostic scans for such interactions are analytically challenging. Furthermore, given the nonlinear effect of genetic load on risk, explicit and strong pairwise interactions may not be the correct hypothesis—gene-based or pathway-based interactions may be more important. Because of limitations in the data, gene-environment interactions were not computed and this area needs study.

The individual roles of DR3 and DR15 haplotypes in SLE risk are well-established. However, in all three ancestries, having two different risk alleles yielded higher SLE risk than having two copies of the same risk allele. This is similar to type 1 diabetes, where heterozygotes for type 1 diabetes-associated haplotypes, DR3 and DR4, have shown higher risk of disease. It is hypothesized that this effect is driven via formation of DQA1 and DQB1 trans-heterodimers. In contrast, SLE risk alleles in DR3 and DR15 stem from divergent ancient haplotypes18; likewise, trans-pairing has not been shown between DQA and DQB in these two haplotypes34,35.

Due to the highly polymorphic nature of HLA alleles and their protein products, it is important to consider high-order relationships among amino acids in three-dimensional space36. Standard regression techniques using amino acids in isolation can be problematic and inappropriate for inference37. To account for higher-order relationships among amino acids, we (1) clustered alleles by protein sequence similarity, (2) compared associations within and between clusters and (3) identified, when possible, amino acids that uniquely distinguished the risk alleles. This approach identified several examples of specific amino acids differentiating risk and protective HLA alleles. For example, the DRB15*01 amino acids −1, 47 and 71 were unique to risk alleles. The combination of Ala71 and Phe47 create a hydrophobic space in the protein binding pocket compared to the alternatives observed (Glu71 and Tyr47; or Arg71 and Tyr47). In addition to antigen binding, there is a vast array of HLA allele-specific properties, including surface expression stability35, influence of DNA methylation38 and DR-DQ heterodimers39. Such findings may help prioritize functional experiments, as we work towards understanding the HLA mechanisms of SLE.

Two major limitations of this study are the comparably fewer non-EA SLE cases and appropriate controls, and the strong EA bias in the Immunochip content. Power calculations using allele frequencies and ORs from EA, and the number of AA cases and controls, yielded 445.5 expected Tier 1 and 2 SNP associations; however, only 64 were observed. Although differences in LD contribute to this result, the highly reduced number of detected associations relative to expected, plus the genetic load analyses, strongly suggest that ancestry-specific and -independent loci contribute to SLE risk. It is imperative to recruit more non-EA populations for genetic studies.

In conclusion, SLE has a strong genetic contribution to risk with ancestry-dependent and ancestry-independent contributions. SLE risk has shared and independent genetic contributions relative to other autoimmune diseases. This genetic risk manifests itself as a nonlinear function of the cumulative risk allele load, a pattern potentially shared across autoimmune and non-autoimmune diseases.

Methods

Study cohort

Multiple studies provided de-identified DNA samples with approval from their respective institutional review boards or ethics committees. These ethics review committees included: Cedars-Sinai Medical Center Institutional Review Board; Central Ethic Committee of Denmark; Centrala etikprövningsnämnden; Comité de Etica de la Investigación de Centro Hospital Universitario Virgen Macarena; Centro de Estudios Reumatológicos. Santiago de Chile; Centro Hospitalar Universitário do Porto, Unidade de Imunologia Clinica e Comissão de Ética; CEPI (Comite de Etica de Protocolos de Investigacion) Institution: Hospital Italiano de Buenos Aires; Cincinnati Children’s Hospital Medical Center Institutional Review Board; Clinical Research Unit, Padua University-Hospital, and Ethics Committee, Province of Padua; Comitato Etico Interaziendale AOU Maggiore della Carità Ethics Committee, Novara, Italy; Comite de Bioetica del Consejo Superior de Investigaciones Científicas; Comité de Docencia e Investigación, Hospital Escuela Eva Perón, Gro Baigorria, Santa Fe, Argentina; Comité de Docencia e Investigación, Sanatorio Parque SA; Comite de etica de la investigacion del HIGA San Martín de La Plata, Argentina; Comité de Ética en Investigación Instituto Nacional de Ciencias Médicas y Nutrición Salvador Zubirán; Comité de Ética en Investigación, Instituto Nacional de Medicina Genómica, Mexico; Comité de Ética en Investigación; Comité de Investigación de la Facultad de Medicina de la UANL y Hospital Universitario ‘Dr José Eleuterio González’; Comite Docencia e Investigacion H.I.G.A. Dr Oscar Alende Mar del Plata; Comitè Ètic d’Investigació Clínica de l’Hospital Clínic de Barcelona; Comités de Ética, Bio Ética y de Investigación. Hospital G. Almenara, Esalud, Lima, Perú; Comites de Ética, Bioetica y de Investigación Hospital Nacional Guillermo Almenara Irigoyen, Lima-Perú; Commission d'Ethique Hospitalo-Facultaire de l'Université catholique de Louvain; Duke University Health System Institutional Review Board; Ethics and Research Committee of Hospital General De Occidente; Fundacion Docencia e Investigacion Hopsital Italiano de Cordoba; Institution of Public Health and Clinical Medicine, Rheumatology, Umeå University, Umeå, Sweden; Institutional Review Board of the University of Puerto Rico Medical Sciences Campus; Institutional Review Board Office Northwestern University; Johns Hopkins University School of Medicine Institutional Review Board; London Central Research Ethics Committee Study sponsor: King’s College London; Medical Ethical Committee (METc) of the University Medical Center Groningen; Medical University of South Carolina Institutional Review Board for Human Research; Northwell Health Human Research Protection Program; Oklahoma Medical Research Foundation Institutional Review Board; omisión Nacional de Investigación Científica y Comisión de Ética en Investigación en Salud, Instituto Mexicano del Seguro Social, México; Regional Ethical Review Board at Karolinska Institutet, Stockholm, Sweden; Regional Ethics Review Board in Linköping; Regional Human Medical Research Ethics Committee of the University of Szeged; SickKids REB; The Institution Review Boards for human research at UCLA; The Local Ethics Committee of the Karolinska University Hospital/Karolinska Institutet, Stockholm Sweden; The University Health Network, Research Ethics Board; Institutional Review Board for Human Use University of Alabama at Birmingham; UC Davis Institutional Review Board; UCSF Human Research Protection Program Institutional Review Board; UHN REB; University Health Network Research Ethics Board and by the local ethics boards of the CaNIOS investigators at the following centres: Montreal General Hospital, St Josephs’ Heath Centre, Winnipeg Health Science Center, Queen Elizabeth II Health Sciences Centre, Ottawa Hospital, Hopital Notre-Dame, Calgary Health Sciences Centre, Centre Hospitalier Universitaire de Sherbrooke, and Hopital Maisonneuve-Rosemount; University Hospital of Gran Canaria Doctor Negrin Research Ethic Committee; University of Chicago Institutional Review Board; University of Southern California Health Sciences Institutional Review Board; University of Texas Southwestern Medical Center Institutional Review Board; Uppsala Ethical Review Board; Wake Forest University School of Medicine Institutional Review Board. All study participants provided written consent prior to study enrolment at the institution where the samples were collected. All SLE cases in this study were required to meet at least four of the eleven American College of Rheumatology classification criteria for SLE40,41.

Genotyping and quality controls

Samples were genotyped on the custom-designed Immunochip Illumina Infinium Assay9 according to Illumina’s protocols, using the Illumina iScan scanner at the following centres: Oklahoma Medical Research Foundation, University of Texas Southwestern, HudsonAlpha Institute for Biotechnology, North Shore-LIJ Health System’s Feinstein Institute for Medical Research. Intensity data were generated for all samples and sent to the Oklahoma Medical Research Foundation for genotype calling using OptiCall42. OptiCall default options were used with one exception: the ‘-nointcutoff’ option was included to allow removal of intensity outliers. Subsequent genotype clusters were viewed against their intensity data using Evoker43. Genotype calling was completed in four batches, keeping samples genotyped at the same center in the same batch. Batches were designed to include samples of multiple ancestries when possible to improve rare variant calling. The ancestry breakdown for the batches was: Batch I was 15% European ancestry (EA), 7% African American ancestry (AA), 55% Asian ancestry (ASA), 23% Hispanic ancestry (HA); Batch II was 44% EA, 18% AA, 1.4% ASA, 36% HA; Batch III was 48% EA, 38% AA, 1% ASA, 13% HA; and Batch IV was 92% EA, 8% AA. Some samples called with the SLE Immunochip study samples were used for other Immunochip studies.

Samples were excluded if their call rates were <98% across SNPs that passed quality control filters. Duplicates and first-degree relatives were removed, retaining the sample with the highest call rate. The Immunochip does not have sufficient markers in the non-pseudoautosomal regions of chromosome X to reliably complete gender checks. Admixture estimates were computed using the program ADMIXTURE44. HapMap phase 2 individuals (CEU: Utah residents with ancestry from northern and western Europe; YRI: Yoruba in Ibadan, Nigeria; CHB: Han Chinese in Beijing, China) as anchoring populations. To facilitate testing for association between rare variants and SLE, and to improve multilocus modelling in regions of linkage disequilibrium (LD) among SNPs, a factor analysis was computed on the admixture estimates using principal component extraction and varimax rotation45. The resulting factors are orthogonal (independent) and thereby remove collinearity among the admixture estimates when used as covariates in linear models. Reduced collinearity should facilitate more robust analysis of rare variants. In addition, principal component (PC) analysis was computed using Eigensoft v4.2 (refs 46, 47) including HapMap phase 2 individuals (CEU, YRI and CHB) as reference populations. Both the admixture and PC analyses were completed using a subset of SNPs generated by removing SNPs in LD (r2>0.2), with minor allele frequency (MAF)<0.01, or with low call rate (<95%).

The admixture estimates and PCs were used to identify and remove genetic outliers. A SNP was removed from the primary analysis if it had an overall call rate <95%, exhibited significant differential missingness between cases and controls (P<0.05), had significant departure from Hardy-Weinberg equilibrium expectations (P<1 × 10−6 in cases, P<0.01 in controls) or a cluster separation score <0.40. SNPs violating the above Hardy-Weinberg equilibrium thresholds were retained if there was convincing evidence of association at SNPs in linkage disequilibrium (LD) and the cluster plots indicated that the pattern was not due to poor genotype calling. Primary inference was based on SNPs with MAF ≥0.01. Finally, >10,000 SNP cluster plots were visually examined, including all SNPs reported, to remove results potentially based on poor genotyping.

To provide an estimate of the number of independent tests for multiple comparisons adjustment, the SNPs were LD pruned, r2<0.20, within each ancestry. The union of these SNPs across ancestries was 46,744 uncorrelated SNPs, yielding a Bonferroni threshold of P<1.06 × 10−6.

Statistical analysis

Regions in figures and tables are named by the genes bounding the regions of association or regions of significance for other statistical test, unless the literature strongly implicated a specific gene.

To test for an association between a SNP and case/control status within an ancestry, a logistic regression analysis was computed adjusting for admixture factors as covariates. Primary inference was based on the additive genetic model unless there was significant evidence of a lack-of-fit to the additive model (P<0.05). If there was evidence of a departure from an additive model, then inference was based on the most significant of the dominant, additive, and recessive genetic models. The additive and recessive models were computed only if there were at least 10 and 30 individuals homozygous for the minor allele, respectively. These tests of association were computed using the SNPGWA version 4.0 module of SNPLASH (https://www.phs.wakehealth.edu/public/bios/gene/downloads.cfm). For ancestry-specific analysis of the X chromosome, the data were first stratified by gender and then meta-analysed using the weighted inverse normal method (weighted by sample size). The genomic control inflation factor (λGC) was calculated using a set of SNPs included on the Immunochip for a study investigating the genetic basis for reading and writing ability. The resulting λGC was scaled to 1,000 cases and 1,000 controls to standardize comparisons across populations and studies.

Three tiers of statistical significance are reported. Tier 1 includes those SNPs that meet the literature-motivated genome-wide threshold of 5 × 10−8. Tier 2 includes those SNPs that are not Tier 1 SNPs, but have a P value for association less than 1 × 10−6. Tier 3 includes those SNPs that do not meet criteria for Tiers 1 or 2, but meet a genome-wide Benjamini–Hochberg false discovery rate48 adjusted P value threshold of 0.05. The Tier 2 threshold meets the strict Bonferroni criteria for the number of uncorrelated SNPs (r2<0.20).

Ancestry-specific logistic regression models were computed to test for evidence of interactions among all pairs of SNPs that had BH-FDR adjusted P value <0.05. Each logistic model contained the admixture factors, the two SNPs, and their centred cross-product term, with the latter term tested using the likelihood ratio test implemented in the Intertwolog module in SNPLASH. To adjust for the number of interactions tested, Bonferroni and BH-FDR adjusted P values were computed. To test for ancestry-specific gender-by-SNP interactions, a case-only autosomal scan was computed; here, gender was the outcome and admixture factors and SNP were the predictors. To adjust for the number of tests computed, the BH-FDR adjusted P values from the likelihood ratio test were computed for each SNP that passed quality control.

To determine how many distinct associations were within a genomic region, a manual stepwise procedure (that is, forward selection with backward elimination, entry and exit criteria of P<0.001) was computed.

For the transancestral meta-analyses, three ancestries were examined for association and meta-analysed to better isolate shared SLE-risk loci by leveraging their LD pattern differences. For each SNP, a nonparametric meta-analysis, weighted inverse normal method (weighted by sample size), was computed as implemented in METAL49. Regions of association were visually examined and tests of heterogeneity of the odds ratio were computed. Thus, for each region, ancestry-specific and meta-analytic tests of association and tests of heterogeneity are reported. The transancestral patterns of association and LD were visualized using LocusZoom50. Results from the weighted inverse normal method were compared to random effects meta-analyses and results of the regions were comparable.

Classical HLA alleles at HLA-A,-B,-C,-DPB1,-DQA1,-DQB1 and -DRB1 were imputed using the program HIBAG51. HIBAG uses an ensemble classifier and bagging technique to arrive at an average posterior probability. Unlike alternative imputation software such as BEAGLE52, HLA*IMP53 and SNP2HLA54, HIBAG did not require training data for any of our three cohorts, as it provides multiple ancestry reference panels (European, African, Hispanic and Asian). This, combined with its accuracy rates being comparable to other approaches51, made HIBAG an ideal method for HLA imputation in our EA, AA, and HA cohorts. To account for imputation uncertainty, the allele dosage was utilized for all analyses. To filter out the lowest frequency alleles, a minimum best guess allele count of 10 was required in either the cases or controls for each allele, in each cohort.

For analysis of classical HLA alleles, single-allele associations were evaluated using logistic regression under the additive model and accounting for imputation uncertainty via allelic dose. To account for population substructure, cohort-specific factors were used as covariates (EA: factors 1–4; AA: factors 1–3; HA: factors 1–2) in each analysis. Meta-analysis was completed for any allele that had a single-allele analysis in at least two cohorts. Evidence of association from each cohort was combined using the weighted inverse normal method via METAL49 and tests for heterogeneity of the odds ratio were computed.

To build multi-locus ancestry-specific models of classical HLA alleles for case/control status of SLE, stepwise regression models were computed. Stepwise logistic modelling (forward selection with backward elimination) was computed using all of the classical HLA alleles that met the QC criteria, including requiring at least a count of 10 alleles from the best guess allele count cross the individuals within an ancestry. The entry and exit criteria were set to P<0.01 for each of the three cohorts. As in the single-allele analysis, the logistic models tested for an additive effect of the alleles and accounted for imputation uncertainty via allelic dose.

To evaluate and compare classical HLA allele associations across the three cohorts, the results from the single-allele and multilocus modelling were visualized in the context of classical HLA protein sequence similarity. Protein sequences for all observed HLA-imputed alleles were retrieved from the EMBL-EBI Immunogenetics HLA Database55. Sequences within an HLA-gene were aligned using ClustalOmega56. Unrooted phylogenetic trees for each of the HLA loci were then generated by Clustal-W2 via the aligned amino acid sequences. The neighbour-joining method, a distance matrix method, utilized a Markov chain of nucleotide or amino acid substitution57. The neighbour-joining method uses this distance information to iteratively evaluate all pairings of neighbours in order to construct a tree that minimizes the branch length at each stage of clustering58. The resulting trees were visualized using Dendroscope59. All results from the single-allele and multilocus classical HLA associations from the three cohorts were graphically displayed on the unrooted trees.

A second set of ancestry-specific single-SNP analyses was computed across the HLA locus and surrounding region, while adjusting for the primary SLE-associated HLA risk alleles from the stepwise modelling. The logistic regression model was computed, as above, considering the fit to the three genetic models (dominant, additive, recessive); the additive model required at least 10 homozygotes for the minor allele, while the recessive model required at least 30. The meta-analysis of these results was computed using METAL.

The Wald tests for HLA-by-SNP and HLA-by-gender interactions were computed using logistic regression models that adjusted for admixture factors and included both the main effects of the HLA allele and SNP (or gender) and their centred cross product as the multiplicative interaction term.

To test whether there was a difference in SLE risk between individuals homozygous for the same risk allele versus heterozygous for two different risk alleles, a Wald test from a logistic regression model was computed adjusting for admixture.

To examine ancestry of associated SLE risk alleles, genotyped SNPs from the population-specific (Tier1 and Tier 2) and the meta-analysis (primary and secondary) tables were compiled into a list of 205 unique SNPs. For evaluation, only SNPs of good quality across the three cohorts were retained. These criteria left 181 SNPs for comparison. In cases, admixture proportions of CEU and YRI were calculated using ADMIXTURE and then the average proportions were tallied for each cohort. Within each of the three populations and for each SNP, the risk allele's average admixture was computed. The resulting risk allele average admixture proportion was compared to the overall average sample admixture proportion in cases by computing the difference between risk allele and sample admixture proportion averages.

To evaluate the SLE-risk allele genetic load, the EA samples were partitioned into two groups: training (the entire EA sample minus 2,000 cases and 2,000 controls randomly chosen from the full EA cohort) and testing (the aforementioned 2,000 cases and 2,000 controls). In the training samples, the single SNP association and stepwise analyses were repeated to obtain a training set of SNPs that had BH-FDR adjusted P-value <0.05. From these results, the EA SLE-risk genetic load was calculated for each individual as the count of risk alleles from the training SNPs. Specifically, we define the EA SLE-risk allele genetic load as:

where, GRSi is the genetic risk score for individual i; γk is the beta coefficient for the kth SNP association with SLE and serves as the weight for that risk allele; RAk is the number of risk alleles for the kth SNP (0, 1, 2); and N is the number of SNPs. By definition of parameterizing relative to the risk allele, γk>0 for all k. The EA SLE-risk genetic load was computed for AA, HA, and the EA testing samples. Individuals whose genetic load (risk allele count) was in the lower 10% of the count distribution were used as the reference sample. A logistic regression model, including admixture factors as covariates, computed the odds ratio comparing the reference sample to samples within a moving window of 20 unweighted risk allele counts for the unweighted analysis and moving window of 4 for the weighted analysis). For example, a logistic model compared the risk of SLE for those in the lowest 10% to those whose risk allele counts ranged from 940 to 960 in the unweighted analysis. The next model and odds ratios were then computed, sliding the allele count up one (for example, 941–961). A plot of these odds ratios for moving windows of 20 counts was constructed to illustrate the pattern. The corresponding plot of the log(OR)=β from the genetic load association with SLE was generated to show that the nonlinearity was not due to the scale; that is, it documents a departure from linearity on the logit scale. A similar approach was completed for a weighted risk allele count, where each risk allele was weighted by the natural logarithm of the odds ratio from the EA SNP association analysis. Plots of the odds ratio effect of the EA genetic load (weighted and unweighted) were generated for AA, HA and the independent EA set.

Finally, for each ancestry an admixture-adjusted regression model was computed to test whether genetic load was associated with age of SLE onset. For ease of interpretation, the strength of the association was reported as the Spearman’s rank correlation coefficient, but the P value is from the admixture-adjusted linear regression model.

Functional annotation analysis

To identify eQTLs for SLE-associated SNPs, all 1,000 Genomes SNPs in LD with the SLE-associated SNP were identified using SNAP60. Specifically, LD was computed using the CEU (for EA and HA) or YRI (for AA) data with an r2≥0.5 for Tier 1 and 2 SNPs. SNPs and their proxies were then queried in a data set downloaded from the eQTL Browser (http://eqtl.uchicago.edu/cgi-bin/gbrowse/eqtl/; Pritchard lab, University of Chicago) and the GTEx Portal (http://www.gtexportal.org). The eQTL Browser contains eQTL data surveyed from 17 eQTL studies, and the Blood eQTL Browser61. The GTEx Portal is a comprehensive resource, with eQTL data from 44 different tissues. When multiple proxies existed for the same eQTL (that is, same SNP and same gene), only the proxy with the lowest P value was retained.

RegulomeDB is a database that annotates SNPs with known and predicted regulatory elements (eQTLs, DNAase hypersensitivity, binding sites of transcription factors) in the intergenic regions of the human genome22. It includes high-throughput, experimental data sets from GEO, the ENCODE project, published literature, as well as computational predictions and manual annotations to identify putative regulatory potential and identify functional variants22. The variants associated with SLE (identified in Tier 1 and 2 in any ancestry cohort) were queried in RegulomeDB.

HaploReg v2 is a tool for exploring annotations of the noncoding genome at variants on haplotype blocks23 and uses LD information from the 1,000 Genomes Project Phase 1 individuals. It analyzes sets of SNPs for an enrichment of cell type-specific enhancers, and includes all dbSNP build 137 SNPs, predicted chromatin state in nine cell types, conservation across mammals, motif instances from ENCODE experiments, enhancer annotations on 90 cell types from the Roadmap Epigenome Mapping Consortium and eQTLs from the GTEx eQTL browser23. The query was performed using default settings, including LD calculations based on the 1,000 Genomes Phase 1 EUR individuals, and epigenome data from both the ENCODE and Roadmap Epigenome Mapping Consortium projects.

SNPs associated with SLE (Tiers 1 and 2) were annotated with the eQTL data and HaploReg v2 (ref. 23) to prioritize those with the highest biological potential. The top summary gene scores were summed across individual criteria (presence of an eQTL, presence of a nonsense or missense variant, promoter and enhancer status in a lymphoblastoid B-cell line (B-LCL), the presence of a DNase hypersensitivity site in any of five immune-related cell lines, presence of a conserved region, the presence of any bound protein, and transcription start site and enhancer status in any of 15 immune cell types), in the haplotype block of each SNP. In the calculation of the biological scores, each functional annotation was given a weight according to their regulatory potential. A score of ‘3’ was given to SNPs in an LD block with any variant that mapped within an active or poised TSS in any of 15 immune cell types, was an eQTL, was non/missense, or mapped within an active promoter in a B-LCL. A score of ‘2’ was given to SNPs in an LD block with any variant that mapped within an active upstream flanking TSS in any of 15 immune cell types or mapped within a conserved region. A score of ‘1’ was given to SNPs in an LD block with any variant that mapped within a weak TSS or any enhancer in any of 15 immune cell types, mapped within a weak promoter or weak enhancer in a B-LCL, mapped within a DNase hypersensitivity site in any of 5 cell lines, or had any bound protein. The sum of these annotations resulted in a final biological score, ranging from zero to fifteen.

For each of the 146,111 (145,278 unique) SNPs that met quality control standards in at least one population, the flanking base pairs were identified using the UCSC reference genome (build 37). Once strand alignment was confirmed between the Immunochip and UCSC reference genome, it was evaluated whether either (or both) of a SNP’s alleles created a CpG site in the 5′-3′ direction.

Data availability

The summary data are available at www.immunobase.org. Individual genotype data, consistent with the respective Institutional Review Board approval and subject consent, are available from the corresponding authors.

Additional information

How to cite this article: Langefeld, C. D. et al. Transancestral mapping and genetic load in systemic lupus erythematosus. Nat. Commun. 8, 16021 doi: 10.1038/ncomms16021 (2017).

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.