Common variants at PVT1, ATG13-AMBRA1, AHI1 and CLEC16A are associated with selective IgA deficiency.

Selective immunoglobulin A deficiency (IgAD) is the most common primary immunodeficiency in Europeans. Our genome-wide association study (GWAS) meta-analysis of 1,635 patients with IgAD and 4,852 controls identified four new significant (P < 5 × 10-8) loci and association with a rare IFIH1 variant (p.Ile923Val). Peak new variants (PVT1, P = 4.3 × 10-11; ATG13-AMBRA1, P = 6.7 × 10-10; AHI1, P = 8.4 × 10-10; CLEC16A, P = 1.4 × 10-9) overlapped with autoimmune markers (3/4) and correlated with 21 putative regulatory variants, including expression quantitative trait loci (eQTLs) for AHI1 and DEXI and DNase hypersensitivity sites in FOXP3+ regulatory T cells. Pathway analysis of the meta-analysis results showed striking association with the KEGG pathway for IgA production (pathway P < 0.0001), with 22 of the 30 annotated pathway genes containing at least one variant with P ≤ 0.05 in the IgAD meta-analysis. These data suggest that a complex network of genetic effects, including genes known to influence the biology of IgA production, contributes to IgAD.

IgAD is the most common primary immunodeficiency and is defined by serum IgA level <0.07 g/L.1 The prevalence in Europeans is 1:600.1 Secretory IgA (sIgA) is important for mucosal immunity and gut commensalism2,3, and clinical features of IgAD include recurrent mucosal infections. In IgAD, B cells fail to terminally differentiate into IgA+ plasma cells, however IL-21 can restore B cell IgA production in vitro4,5. IgAD is strongly associated with HLA6 and aggregates in families with autoimmunity.7 The prevalence of Celiac disease is 35 times higher in IgAD patients whereas the prevalence of systemic lupus erythematosus (SLE) and type 1 diabetes (T1D) is 10 times higher.8

Ferreira et al.6 previously imputed HLA-B, HLA-DRB1 and HLA-DQB1 in IgAD cases and controls (>2,700 individuals).9 The primary signal mapped to HLA-DQB1*02 (OR = 2.8, P = 7.7x10−57), due to combined independent effects of the HLA-B*08:01-DRB1*03:01-DQB1*02 and HLA-DRB1*07:01-DQB1*02 haplotypes.9 There was a secondary signal at HLA-DRB1*01:02 (OR = 4.28, P = 5.86x10−17) and a protective effect for HLA-DRB1*15:01 (OR = 0.13, P = 2.24x10−35).9 HLA-DQB1*02:01 is protective for IgA nephropathy (OR = 0.71, P = 2.61x10−13). 10 None of the non-MHC IgAD loci overlapped with IgA nephropathy.10

A previous GWAS identified the common allele IFIH1 Thr946Ala (OR = 0.62, control allele frequency = 0.39) as the first and, to date, only non-HLA genome-wide significant IgAD locus.9 Thr946Ala is also protective for T1D11, SLE12, psoriasis13 and vitiligo (see URLs). IFIH1 encodes MDA5, a cytosolic receptor that recognizes dsRNA and initiates interferon pathway activation.

To expand our understanding of IgAD risk, we studied four new IgAD cohorts, and performed a GWAS meta-analysis of ~9.5M SNPs in 1,635 cases and 4,852 controls (Table 1). Genotypes for untyped markers were imputed for each cohort separately (1000 Genomes Project) and genotypes for variants fully typed in the entire cohort. Up to four controls per new case were iteratively selected based on ancestry eigenvectors14 to minimize population substructure (Online Methods). Genomic inflation factor values were minimal (Table 1), indicating that population substructure was adequately addressed. Association analyses were conducted with logistic regression (additive model), accounting for genotype uncertainty and using ancestry eigenvectors as covariates (Online Methods).

Table 1

Case-Control Cohorts for IgAD GWAS.

Cohortsa	Cases (N)	Controls (N)	Total (N)	Significant Ancestry Eigenvectors (N)	Genomic Inflation Factor (λGC)b	SNPs Genotyped (N)	SNPs after Imputation (N)	Median Info Score
Meta	1,635	4,852	6,487			556,344a	9,464,381
2016 (New)
Swedish	483	1,932	2,415	2	1.03	423,694	9,464,381	0.988
Spanish	150	230	380	0	1.04	339,552	8,534,763	0.974
Italian	91	364	455	4	1.03	218,770	8,365,602	0.948
Czech	151	602	753	4	1.05	112,822	6,282,267	0.929
2010
Swedish	421	1,080	1,501	5	1.03	289,843	8,765,152	0.979
Finnish	86	344	430	3	1.01	314,756	8,723,593	0.984
Spanish	253	300	553	3	1.00	538,800	8,794,151	0.985

Variants genotyped in >1 cohort.

The genomic inflation factor listed is calculated from genotyped variants. However, genomic inflation factors for each cohort were also estimated for imputed SNPs separately, and results did not differ (e.g. maximum genomic inflation factor for imputed variants was 1.06).

The strongest association was with the MHC region, 2.7kb upstream of HLA-DQA1 (P = 3.3x10−92) (Table 2). Protective association between IgAD and IFIH1 Thr946Ala9 was confirmed (P = 3.7x10−15), and association with the rare loss of function IFIH1 Ile923Val variant (P = 2.6x10−8) was subsequently identified (Supplementary Table 1). Val923 has previously been shown to be protective for T1D11 and psoriasis13 and abrogates interferon signaling.15 A rare gain of function IFIH1 Arg779His variant has also been reported in an IgAD case with SLE and a type I interferon signature.16

Table 2

Genome-wide Significant Results (P < 5x10−8) for IgAD GWAS.

Closest Gene(s)	Varianta	Position (hg19)	Minor Allele	% Frequency in Swedish Sample		P	FDR q values	OR	Type	Immune Diseases Sharing this Locusb
				Cases	Controls
HLA-DQA1 c	rs116041786	6:32602396	C	14.4	38.7	3.3e–92	1.0e–87	0.38	Intergenic	Celiac Graves’ IBD Sjögren’sd
IFIH1	rs1990760	2:163124051	G	29.3	37.3	3.7e–15	2.7e–12	0.70	Missense	T1D SLE Psoriasis Vitiligoe
PVT1	rs11299600	8:129204573	1 bp deletion	18.9	25.3	4.3e–11	2.4e–8	0.73	Intergenic	RAf
ATG13 – AMBRA1	rs4565870	11:46349869	C	28.9	23.6	6.7e–10	3.5e–7	1.38	Intergenic
AHI1	rs7773987	6:135707486	C	51.3	44.6	8.4e–10	4.3e–7	1.30	Intron	MS
CLEC16A	rs34069391	16:11161214	1 bp insertion	14.3	19.6	1.4e–09	6.9e–7	0.71	Intron	T1D MS PBCg

Peak variants were imputed, except for IFIH1 and CLEC16A, which were genotyped in all but the Czech cohort. Three loci had suggestive association in our previous GWAS: PVT1 (P = 4x10−6), DGKZ (P = 2x10−6) and CLEC16A (P = 2x10−7).

Autoimmune diseases reporting a genome-wide significant variant (P < 5x10−8) in Immunobase and the GWAS Catalog (see URLs) with at least a modest effect (OR ≥ 1.1) in the same direction as IgAD results and on the same haplotype as a peak IgAD variant (r > 0.45).12 IBD, inflammatory bowel disease; MS, multiple sclerosis; SLE, systemic lupus erythematosus; PBC, primary biliary cirrhosis; RA, rheumatoid arthritis; T1D, type 1 diabetes.

Genetic association of IgAD with HLA class II is well-established, particularly with HLA-DQB1*02:01 and HLA-DRB1*01:02 and the highly protective HLA-DRB1*15:01 allele.6 Few of the cases (2.1%) were carrying two copies of rs116041786*C (versus 15.7% of controls).

MS, RA, SLE, T1D, vitiligo, allergy and asthma also have some associations reported at this locus.

IBD is also associated with this locus.

MS, celiac disease, eczema and allergy also have some associations reported at this locus.

Celiac disease and asthma also have some associations reported at this locus.

Our GWAS identified four new significant (P < 5x10−8) loci: PVT1, ATG13-AMBRA1, AHI1 and CLEC16A (Figure 1, Table 2). Peak novel IgAD variants were (1) a protective one bp deletion 91kb downstream of PVT1 (P = 4.3x10−11); (2) a variant 289kb upstream of ATG13 (P = 6.7x10−10); (3) an intronic AHI1 variant (P = 8.4x10−10); and (4) a protective one bp intronic insertion in CLEC16A (P = 1.4x10−9) (Table 2, Supplementary Table 2, Supplementary Figures 1-6). 28 additional loci contained at least one variant with 5x10−8 < P < 1x10−5, and half (denoted with *) were known autoimmune loci: FAS*, GATA3*, MYO9B, CD86*, BCL6, IKZF2, PTGER4*, IKZF3*, EGR2*, FCRL3*, IL2RA*, PTPN2*, IL2*/IL21*, CDH23*, MAST4, FOXP1, MIR605, LINC01098, DRMT1, LINC00299, C22orf42, FAM171A1, BACH1, ZCCHC24, NCKAP5, LOC100505887, TMEM72 and CD28* (Supplementary Table 3, Supplementary Figure 7).

Figure 1

Genome-Wide Significant Loci in the IgAD Meta-Analysis.

Manhattan plot of P-values from a genome-wide meta-analysis of ~9.5M variants in 1,635 IgAD patients and 4,852 controls. The vertical line indicates genome-wide significance (P < 5x10−8). Loci with genome-wide significant results are labeled; red labels indicate novel genome-wide significant loci.

The heritability of IGAD in the Swedish cohort based on genome-wide imputed variants was 0.49 (standard error [S.E.] = 0.047). This estimate was reduced from 0.49 to 0.39 (S.E. = 0.059) when we conditioned on the five peak non-MHC variants, the peak MHC variant (rs116041786), and the peak MHC variant after conditioning on the peak MHC variant (rs116350876). When we excluded the MHC region altogether from the input, the estimated heritability dropped to 0.14 (S.E. = 0.054).

As one approach to compile a list of potential causative gene(s) and allele(s) in the new loci,17 we cross-referenced the peak variants plus 160 correlated variants (r ≥ 0.7) (44 variants in PVT1; 21 variants in ATG13–AMBRA1; 65 variants in AHI1; and 30 variants in CLEC16A; total N = 164), against epigenetic data, eQTLs, DNase I hypersensitive sites, transcription factor binding sites (ChIP-Seq), consensus motifs, GWAS loci, and promoter regions using RegulomeDB18, ENCODE19, Roadmap Epigenomics, GEO, the WashU EpiGenome Browser, ImmunoBase, and the GWAS Catalog (see URLs). The RegulomeDB catalogues lists regulatory evidence for ~60M variants, and about 5% (~3M) are categorized as likely to affect binding (score = 2) or likely to affect binding and linked to expression of a gene target (score = 1).

In our dataset, 13.1% (21/160) of the correlated variants were likely to affect binding (RegulomeDB score < 3), and were located in regions of open chromatin in an immune cell (as determined by DNase hypersensitivity) and/or in transcription factor binding sites (as determined by DNase footprints, ChIP-Seq, and/or binding site motifs) (Supplementary Table 4). Many of these variants (15/21) were also located in histone marks for active transcription (H3k4me3, H3k27ac or H3k4me1) in T regulatory cells (Treg) and all of the novel loci, except for AHI1, encompassed at least one Treg enhancer and active transcription start site. A comparison of the overlap between eQTL data from lymphoblast B cells and IgAD meta P-values using the Bayesian statistical framework Sherlock identified three cis eQTLs in AHI1 and DEXI associated with IgAD (Supplementary Table 4).

The PVT1 locus contains no protein coding genes, and the lncRNA PVT1 appears to be the most likely causative genetic element. PVT1 has been shown to be important for expression and copy number increase of MYC in tumors.20 The peak PVT1 variant was in moderate LD with four potential regulatory variants. The most interesting of these was rs7001706 (r= 0.70), an intronic variant that lies in a FOXP3 transcription factor binding motif and in an H3K4me1 histone mark in Tregs. When we conditioned the regression analyses on either the peak variant or the potential regulatory variant rs7001706, they accounted for part but not all of the locus effect (data not shown), indicating that PVT1 contains more than one independent signal.

AHI1 is highly expressed in hematopoietic cells and overexpressed in blood cancers.21 AHI1 stabilizes BCR-ABL in leukemia cells by recruiting JAK2, and can regulate phosphorylation of JAK2-STAT5.22 JAK inhibitors are useful in treating rheumatoid arthritis (RA), psoriasis and alopecia areata. The peak variant was in LD with 8 eQTLs, all on the same haplotype. Rs2064430 (r= 0.85) had the strongest evidence for regulatory binding (RegulomeDB score = 1d) (Supplementary Table 4). A comparison of the overlap between eQTL data from lymphoblast B cells and IgAD Pmeta using the Bayesian statistical framework Sherlock highlighted two cis eQTLs in AHI1 (rs2179781, rs9647635) with strong positive Log Bayes Factor (LBF) scores (6.83, 6.63) (Supplementary Table 4). Of the variants highlighted in Supplementary Table 4, rs2179781 had the smallest Pmeta value. Rs2179781*A was associated with reduced risk for IgAD and reduced AHI1 expression (Supplementary Figure 8).

The peak SNP upstream of ATG13 (autophagy 13) and AMBRA1 (autophagy/beclin-1 regulator 1) was correlated with two likely regulatory variants (Supplementary Table 4). ATG13 and AMBRA1 physically interact 23 and share the same chromatin interaction domain. Autophagy plays a role in autoimmunity24, plasmablast differentiation25, and immunoglobulin production in plasma cells.26 ATG13 and AMBRA1 are the most likely risk genes in this interval, but the causal variant remains to be determined.

CLEC16A contains several eQTLs for DEXI,27 and T1D-associated CLEC16A variants28 are located in an intron that shows chromatin interaction with the DEXI promoter, as determined by 3C27 and Hi-C37,38 mapping. The peak variant was in LD with 7 likely regulatory variants (Supplementary Table 4), including rs35300161, rs34972832 and the DEXI eQTL rs17806299 (r = 0.80). Rs35300161 and rs34972832 are present in open chromatin in germinal center and naïve B cells (Supplementary Table 4) and in an enhancer in various T cell types and tissues. A comparison of the overlap between eQTL data from lymphoblast B cells and IgAD Pmeta with Sherlock highlighted one cis eQTL in DEXI (rs741175) with a strong positive Log Bayes Factor (LBF) score (6.16) (Supplementary Table 4).

Clec16a knockdown mice show reduced numbers of B cells,29 and are protected from autoimmunity.30 Further work is required to determine whether CLEC16A or DEXI is the relevant gene in this interval.

We next applied a pathway analysis to the meta-analysis results using genomic randomization (PARIS)31 to estimate empirical significance (N = 10,000 randomizations) and an alternate pathway algorithm based on LD-independent intervals for replication (INRICH).32 PARIS reduces bias from gene size, pathway size, LD blocks and SNP coverage. Preliminary testing with PARIS v1.1 using imputed data has not been completed and thus, only the GWAS chip content was selected for this experiment. Using PARIS, we grouped variants into LD features and single SNPs into linkage equilibrium (LE) features. The LD and LE features were grouped by KEGG pathways and also by our user-defined Treg signature pathway. The features in each pathway were permuted 10,000 times with a randomly selected set of features of similar size. The empirical P is based on the number of features with a P < 0.05 in a pathway, compared to the number of features with a P < 0.05 in the permuted pathway. IgAD was associated (empirical P < 0.0001) with 7/221 (3%) KEGG pathways (total genes queried N = 5,701) (Supplementary Table 5).

Of particular interest, one of the top pathways identified was the intestinal immune network for IgA production (empirical P < 0.0001), with 22/30 genes in the pathway containing variants with a P value < 0.05 in the IgAD meta-analysis (Figure 2, Supplementary Table 6). These include genes encoding inflammatory and immune regulatory cytokines (IL-2, IL-4, IL-6, IL-10, IL-15, APRIL, BAFF), key cell surface molecules on B and T cells (TACI, CD40, IL-15R, CD28, ICOS, ICOSLG), molecules essential for the homing of lymphocytes to the gut (ITGA4, CCR9, MADCAM-1, CXCL12, CCL28, LTβR, MAP3KI4) and a DNA deaminase required for B cell class switch recombination (AICDA). These data suggest that the genetic contribution to IgA deficiency includes many genes known to regulate IgA production.

Figure 2

KEGG Pathway for IgA Production Associated with IgAD

A GWAS-based pathway analysis of IgAD meta-analysis P-values was performed of all 221 KEGG database pathways using PARIS and INRICH software (see Methods). A top identified pathway (P < 0.0001 based on 10,000 permutations) was the “Intestinal Immune Network for IgA Production” where 22/30 annotated pathway genes contained variants with a P value < 0.05 in the IgAD meta-analysis. The colors of the boxes around the gene names correspond to the smallest P value observed in the IgAD meta-analysis within that gene: blue, P < 1x10−5; green, 1x10−5< P < 1x10−3; yellow, 1x10−3< P < 0.01; red, 0.01 < P < 0.05, white P ≥ 0.05.

Similar to RA and T1D GWAS SNPs, which are enriched for active histone marks (H3k4me3) in Tregs17 and lymphoid gene enhancers33, respectively, IgAD variants were enriched for active histone marks and enhancers in Tregs. When we applied a set of Treg signature genes (27 autosomal genes that are significantly upregulated in FOXP3+ Tregs versus conventional FOXP3– CD4+ T cells34) as a user-defined pathway in a PARIS pathway analysis, we found that the Treg signature pathway was significantly enriched for association with IgAD (empirical P < 1x10−5) (Supplementary Table 7). FOXP3 is an important transcriptional regulator of Treg differentiation,35 and germline mutations of FOXP3 cause IPEX,36 a disease characterized by impaired Treg function and life-threatening autoimmunity. Tregs in B cell germinal centers of Peyer’s patches have been shown to influence the ability of sIgA to maintain microbial diversity.37

Using GWAS summary statistics for T1D and RA, we tested for shared genetics between these two diseases and IGAD. Applying GPA,38 we found significant shared genetics both between T1D and IGAD, and between RA and IGAD. GPA uses a likelihood ratio test to compare a model of the pairwise summary statistics with no shared genetics to a model with shared genetics. Results indicate significant shared genetics, even when variants in the MHC region were omitted (P < 1.6x10−7).

In conclusion, we identified four new susceptibility loci for IgAD – PVT1, AHI1, AMBRA1–ATG13 and CLEC16A. Interrogation of peak IgAD variants and SNPs in LD revealed 21 putative regulatory variants. Pathway analyses highlighted enrichment for association in the intestinal immune network for IgA production and in Treg signature genes. Further work is needed to validate novel associations via genotyping, to identify causal variants in the new loci, to explore the regulatory role of associated variants through functional studies, and to investigate autoimmune mechanisms that contribute to IgAD pathogenesis.

Online Methods

The GWAS study tested ~9.5M variants in 1,635 IgAD patients (clinical diagnosis defined as serum IgA ≤0.07 g/L with normal IgM and total IgG)39 and 4,852 controls from seven independent case-control cohorts (see Table 1).

Samples

Only IgAD cases without comorbid autoimmune diseases (e.g. Type 1 diabetes, celiac disease) were included in this study. Furthermore, IgAD cases were screened for celiac disease (auto-antibodies to tissue transglutaminase) and excluded from the study if positive.

Genotyping

Illumina, Inc. (San Diego, CA, USA) developed all of the genotyping arrays. We genotyped cases on Omni1-Quad and Omni2.5, Swedish controls (Swedish TwinGene Project)40 on Omni Express, Spanish controls on Omni1-Quad, Italian controls on 1Mb (Milan Hypergenes Study)41 except for 42 Italian controls genotyped on the 550K (New York Cancer Project (NYCP))42, and Czech controls on HumanHap1M, Omni1-Quad, and Omni 2.5 (National Cancer Institute Division of Cancer Epidemiology and Genetics) Imputation Reference Dataset). 64 Czech controls were genotyped on HumanHap500 (NCI CGEMS (Cancer Genetic Markers of Susceptibility) Project). Ferreira et al. genotyped the Swedish cohort43, which included 34 Icelandic cases, on Hap300 (except for 175 cases genotyped on Human610), Finnish cases on 610-Quad, Finnish controls on CNV370 arrays, and the Spanish cohort on 610-Quad v1. See the “SNPs Genotyped (N)” column in Table 1 for the number of SNPs that overlapped between the arrays for each of the 7 cohorts. All SNPs were mapped to build hg19 coordinates using liftOver.

QC Post-Genotyping and Pre-Imputation

Prior to imputation, we removed variants with a genotyping rate <98%, ambiguity (A/T SNPs and C/G SNPs), evidence of deviation from Hardy-Weinberg equilibrium in controls (P < 1x10−4), and MAF < 1x10−6. Subsequently, we removed individuals who were missing >1% of genotypes and who had heterozygosity ±6 standard deviations (s.d.) from the mean. To estimate identity-by-descent (IBD) and conduct principal components analysis (PCA) (EIGENSOFT v3.044), we randomly selected 100,000 LD-pruned variants (LD pruning was based on r > 0.2; excluding insertion/deletions, variants in regions of long range LD45 or in the MHC (chr6:24–36Mb) or in a large inversion on chromosome 8,46 and variants with MAF < 0.05) and included 878 ancestry informative markers.47 We estimated and plotted IBD for all pairs of samples and removed duplicate samples and cryptic related individuals (P(IBD=0)<0.8). For the PCA analysis, we used 5 outlier removal iterations and 10 principal components (PCs) along which to remove outliers during each outlier removal iteration. Any sample exceeding 6 s.d. along one of the PCs was considered a population outlier and was removed. We then reran PCA without the ancestry outliers.

For the four new cohorts, we iteratively selected up to four controls similar to a previously described approach (developed by M.F.S.) to minimize population substructure differences between cases and controls.48 Each control could only be selected once. For the matching we sequentially selected a control that passed QC and best matched the case, considering the EIGENVALUES results from PCs contributing at least 1% of variance determined from PCs with significant P–values (P < 0.05) in the Tracy-Wisdom test44. We estimated percent variance explained by subtracting the eigenvector value of the first non-significant PC from each significant PC, summing this value, and dividing each eigenvector value by the sum. EIGENVALUES for each PC were adjusted for the EIGENVECTOR scale (i.e. adjusted for the relative contribution of each EIGENVECTOR in describing the data). The lowest cumulative sum of the absolute differences (Manhattan distance for each PC after adjustment based on relative EIGENVALUES) between all cases and all controls was scored as the best match.

Subsequently, we reran PCA a third time for each cohort using the final set of case-control samples. PCs contributing at least 1% of variance determined from PCs with significant P-values (P < 0.05) in the Tracy-Wisdom test44 were used as covariates in the association analyses. Significant PCs were also plotted with cases and controls color coded differently for visual inspection. We tested for differential genotyping rates between cases and controls but did not observe any.

Statistical Analysis Before Imputation

We conducted χ tests of association on genotypes for each cohort separately, using only variants that overlapped between the arrays used to genotype in the cohort, and estimated a genomic inflation factor (λGC) 49 based on the median χ for non-MHC genotyped variants adjusted for ancestry (PLINK v1.0749). Prior to meta-analysis and imputation, we plotted P-values in a quantile-quantile (Q-Q) plot (excluding the MHC) and in a Manhattan plot in R.

Imputation

We imputed genotypes for each cohort separately using the 1000 Genomes Project Phased I integrated variant set release as a reference panel (SHAPEIT v2.r72750, IMPUTE2 v2.3.051). We imputed cases and controls from the same cohort together, using only variants that overlapped between the arrays used to genotype the cohort. Genotypes were imputed in 5 MB chunks, excluding centromeres. We used an effective size of the population of 20,000; 80 haplotypes as templates when phasing observed genotypes; and performed 30 MCMC iterations, where the first 10 MCMC iterations were discarded as burn-in.

Statistical Analysis After Imputation

We excluded variants with an imputation info score <0.6 or a MAF<0.01 in the controls or evidence of deviation from HWE (P < 1x10−6), and tested association for each case-control cohort separately with logistic regression (additive model), adjusted for ancestry (SNPTEST v2.5-beta452). The median info score for each cohort was: Swedish 0.988; Spanish 0.974; Italian 0.948; Czech 0.929; Swedish 2010 0.979; Finnish 2010 0.984; and Spanish 2010 0.985. We used a missing data likelihood score test to address genotype uncertainty. We estimated a genomic inflation factor (λGC)49 based on the distribution of P-values for each cohort using the GenABEL estlambda package in R (see URLs). We plotted P-values in a quantile-quantile (Q-Q) plot (excluding the MHC) and in a Manhattan plot in R.

Meta-analysis

We pooled results in a meta-analysis (PLINK v1.0749), and report P values and ORs for a random effects model. We plotted genome-wide results with R v3.0.2 (see URLs) and regional results in LocusZoom v1.3, using LD estimated in PLINK with the 1000 Genomes EUR data (release 20130502 v.5).53 Meta-analysis P values were adjusted for multiple testing with the false discovery rate (FDR) in R.54

Pathway Analysis of IgAD GWAS Results

We conducted a GWAS based pathway analysis on IgAD GWAS meta-analysis P-values after imputation for a subset of variants (~540K variants from the Illumina 550K array) using PARIS (v 1.1.0b)55. PARIS groups SNPs into LD features and single SNPs into linkage equilibrium features. The LD features can extend up to 50KB past gene boundaries. It then groups these features by pathway, and permutes the genomic structure of interrogated pathways to determine the significance of the pathway while accounting for differences in LD, gene size, pathway size and SNP coverage between pathways. The total number of features with a significant P-value (defined as P < 0.05), is compared with the number of significant features in the permuted pathway. The software was tested and validated using genotype data from the Illumina 550K. We only used a subset of variants because PARIS (v 1.1.0b) has been demonstrated to be robust when used on this SNP array (Yaspan, pers. commun.).

Variants with significant P-values (P < 0.05 from the GWAS meta-analysis data) in Kyoto Encyclopedia of Genes and Genomes (KEGG)56 pathways (v. 3/2011) were identified, and the genome was permuted (10K permutations) to estimate pathway significance (P < 0.0001). Out of 10,000 randomizations, none of the randomized pathways had more variants with P < 0.05 than the identified pathways. A follow up pathway-based analysis was conducted with Inrich (v.1.0) using only pathways that were significant in the PARIS analysis (P < 0.0001).57 We identified clumps (LD independent intervals) using P < 0.0001 as the significance threshold for index SNPs and P < 0.05 as the secondary significance threshold for clumped SNPs and r > 0.5 as the LD threshold for clumping (PLINK). There were 72 non-overlapping LD independent genomic intervals of enriched association, based on the subset of meta-analysis P-values used in the PARIS analysis. Pathways that contained at least two non-MHC genes and had both a PARIS P < 0.0001 and an Inrich P < 0.05 were considered statistically significant.

Regulatory Variation

For each of the five non-MHC IgAD loci reaching genome-wide significance (P < 5x10−8), SNPs that were in LD (r > 0.5, EUR) with the peak SNP were identified. The peak SNP was defined as the SNP with the lowest P-value. LD was estimated using whole genome sequence (WGS) data (mean 30x coverage) available for 583 European Americans (PLINK v1.90b20). Significant peak variants, and variants correlated with the peak variants, were queried for known and predicted regulatory DNA elements (RegulomeDB v1.1, dbSNP141).

Overlap between expression QTLs and GWAS

The overlap between IgAD meta-analysis P-values and expression QTL data from lymphoblast B cells 58,59 was assessed using the Bayesian statistical framework Sherlock (see URLs).60

Heritability Estimates

The program GCTA-GREML (see URLs) was used to estimate the variance between the cases and controls by the entire genome, and the variance explained by the IgAD genome-wide significant variants. The heritability of IgAD for the entire genome was estimated with GCTA in the Swedish IgAD cohort (n = 1,214,325 imputed genotypes).61 The other IgAD cohorts were underpowered for heritability estimates (gctaPower). We estimated the heritability of IgAD for the 7 peak IgAD variants by including the 7 peak IgAD variants as covariates. This analysis was then repeated without the MHC region (chr6:24–36 Mb).

Shared Genetics with T1D and RA

The program Genetic analysis incorporating Pleiotropy and Annotation (GPA)62 was used to estimate the extent of the genetic overlap between IgAD and half a million T1D variants63 and between IgAD and 6 million RA variants64. Summary statistics for T1D were publicly available from dbGaP (accession number: phs000180.v1.p1). Summary statistics for RA were from the European subset of a recent meta-analysis of RA64 (see URLs). There were 5,966,608 and 501,001 variants in the IgAD-RA and IgAD-T1D pairwise analysis, respectively. Using these genome-wide summary statistics, we estimated the proportion of shared variants (variants associated to both IgAD and RA, or to both IgAD and T1D). We tested for significance of shared genetics between pairs of diseases using a likelihood ratio test62. GPA estimates the proportion of variants associated to both diseases and uses a likelihood ratio test to assess statistical significance. This analysis was conducted genome-wide, and then repeated with the exclusion of the MHC region (chr6: 24–36 Mb) due to the extensive LD and highly significant variants in this region. After removing variants in the MHC, 5,944,048 and 498,320 variants remained in the IgAD-RA and IgAD-T1D analyses, respectively.

Supplementary Material

Note: Supplementary information is available on the Nature Genetics website.