Multi-ancestry genome-wide association study of 21,000 cases and 95,000 controls identifies new risk loci for atopic dermatitis.

Genetic association studies have identified 21 loci associated with atopic dermatitis risk predominantly in populations of European ancestry. To identify further susceptibility loci for this common, complex skin disease, we performed a meta-analysis of >15 million genetic variants in 21,399 cases and 95,464 controls from populations of European, African, Japanese and Latino ancestry, followed by replication in 32,059 cases and 228,628 controls from 18 studies. We identified ten new risk loci, bringing the total number of known atopic dermatitis risk loci to 31 (with new secondary signals at four of these loci). Notably, the new loci include candidate genes with roles in the regulation of innate host defenses and T cell function, underscoring the important contribution of (auto)immune mechanisms to atopic dermatitis pathogenesis.

Atopic dermatitis (eczema) is a common inflammatory skin disease affecting 15–30% of children and 5-10% of adults[1]. Its pathogenesis involves skin barrier abnormalities and a T-cell-driven cutaneous inflammation. Atopic dermatitis has significant genetic contributions, with heritability estimates of up to 90%[2] in Europeans. The strongest known risk factors are null mutations of the filaggrin (FLG) gene, resulting in epidermal barrier deficiency[3-5]. Genome-wide association (GWA) studies have identified 20 additional loci (10 in Europeans, 8 in Japanese, 2 in Chinese populations), mostly implicated in immune dysregulation[6-12]. Genetic modeling suggests further loci may be identified with well-powered GWAS[13]. We therefore carried out a multi-ethnic meta-analysis of 26 studies comprising 21,399 cases and 95,464 controls imputed to the 1000 Genomes Project Phase 1 reference panel (Supplementary Note 1 & Supplementary Table 1). 15,539,996 variants with ≥1% MAF were analyzed.

A fixed effects meta-analysis of the 22 European studies identified 21 genome-wide significant (p<5×10−8) loci (Table 1, Fig 1, Supplementary Figs 1-4), and a multi-ethnic meta-analysis identified an additional 6 loci with log10 Bayes Factor>6.1, 4 of which (10q21.2, 6p21.33, 11p13, 2p13.3) also showed nominal association in the European analysis (Table 1). These 27 loci included all 11 loci previously associated with atopic dermatitis in Europeans and 5 loci originally reported in Japanese. Three Japanese loci (6p21.33, 10q21.2, 2q12.1) were also strongly associated in the European analysis, whereas two (3q13.2, 11p15.4) may represent Japanese-specific signals (Supplementary Figs 1&2), with the European confidence interval ruling out all but very small effects (OR<1.03, Table 1). Furthermore, a locus originally reported in a Chinese GWAS (20q13.33) showed association in Europeans. We identified 11 novel loci for atopic dermatitis. Four (11q24.3, 10p15.1, 8q21.13, 2p25.1) were previously associated with self-reported allergy[14], and another (8q21.13) with asthma[15]. Two novel variants (5p13.2 and 2p25.1) showed statistically significant evidence of heterogeneity between European and non-European studies (Cochran’s Q p~0.01, Supplementary Table 2). Both showed little evidence for association in non-Europeans (particularly Japanese, Supplementary Fig.2). The CIs also overlapped for all variants when comparing pediatric (defined as onset by age 6) with any-age onset studies (Supplementary Fig.3). Within Europeans there was some evidence of heterogeneity in effect sizes between studies amongst known variants (e.g. 11q13.5 I2=62.9%, p<0.0001; 11p13 I2=55.6%, p=0.0011) but little evidence amogst novel variants (I2 range=0-40%, all p>0.02, Supplementary Fig.2). Nevertheless, studies with phenotype definition based on a dermatological exam tended to report larger effect sizes than studies using self-report (Supplementary Fig.4), which is to be expected, assuming a moderate degree of phenotypic misclassification in the latter. The inclusion of studies utilizing self-report is therefore likely to bias estimates of the effect size towards the null, and this should be borne in mind when interpreting the odds ratios from our study. Given the primary aim of GWA studies is the detection of novel loci, the increase in sample size achieved by including these studies is so large that any potential detrimental effect on statistical power is more than outweighed and the expected direction of bias means there is unlikely to be an issue of spurious findings (corroborated by Supplementary Fig.4).”

Table 1

Discovery Results. The index variant for loci with p<5×10−8 in the European analysis or log10BF>6.1 in the multi-ethnic MANTRA analysis. Previous atopy trait associations with these loci are listed.

				European – fixed effects				All cohorts – MANTRA		Known atopy loci?
Variant	Locus	Nearest Gene‡	EA/OA	N (studies)	EAF	OR (95% CI)	P-value	N (studies)	log10 BF	trait	references
KNOWN LOCI
rs61813875	1q21.3	CRCT1/LCE3E (FLG)§	G/C	93,326 (18)	0.02	1.61 (1.48–1.75)	5.6×10−29	96,419 (20)	25.53	AD	3,4,5
rs10791824	11q13.1	OVOL1	G/A	102,761 (21)	0.57	1.12 (1.09–1.15)	2.1×10−19	116,556 (25)	21.56	AD	9
rs12188917	5q31.1	RAD50/IL13	C/T	102,761 (21)	0.21	1.14 (1.10–1.17)	4.0×10−17	116,554 (25)	17.24	AD,A,IgE	9,18,58
rs6419573	2q12.1	IL18R1/IL18RAP	T/C	102,760 (21)	0.26	1.11 (1.08–1.14)	1.5×10−13	116,557 (25)	18.10	AD,A,AS,SRA	8,14,18,21
rs2212434	11q13.5	C11orf30/LRRC32	T/C	102,761 (21)	0.45	1.09 (1.07–1.12)	4.6×10−13	116,557 (25)	13.02	AD,AS,SRA,AR,A	11,14,15,21,59
rs4809219	20q13.33	RTEL1–TNFRSF6B	C/A	102,760 (21)	0.27	0.90 (0.87–0.93)	7.0×10−13	116,555 (25)	11.98	AD	7,10
rs2918307	19p13.2	ADAMTS10/ACTL9	G/A	100,707 (20)	0.16	1.12 (1.08–1.16)	4.6×10−12	114,504 (24)	12.98	AD	9
rs2041733	16p13.13	CLEC16A	C/T	103,066 (22)	0.55	0.92 (0.90–0.94)	2.5×10−11	116,862 (26)	10.11	AD,A+HF	7,53
rs12730935*	1q21.3	IL6R	A/G	102,760 (21)	0.39	1.08 (1.05–1.11)	6.1×10−11	116,556 (25)	7.15	AD,A	12,15
4:123243592†	4q27	KIAA109 (IL2)§	R/I	102,761 (21)	0.37	1.08 (1.05–1.10)	4.2×10−9	107,119 (24)	7.32	AD,AS,SRA	7,14,21
rs4713555	6p21.32	HLA-DRB1/HLA-DQA1	T/G	91,217 (15)	0.27	0.91 (0.89–0.94)	5.4×10−9	105,014 (19)	10.76	AD,AS,SRA,A	6,8,14,18,21
rs2944542	10q21.2	ZNF365	C/G	102,762 (21)	0.41	0.94 (0.92–0.96)	1.2×10−6	116,559 (25)	7.56	AD	8,10
rs145809981	6p21.33	MICB	T/C	97,697 (19)	0.14	0.91 (0.88–0.95)	1.5×10−6	110,228 (22)	7.33	AD,AS,SRA	8,14,21
rs4312054	11p15.4	OR10A3/NLRP10	G/T	102,760 (21)	0.41	1.00 (0.97–1.02)	0.744	116,556 (25)	7.00	AD	8
rs1249910	3q13.2	CCDC80/CD200R1L	A/G	99,164 (20)	0.34	0.98 (0.96–1.01)	0.137	112,960 (24)	6.86	AD	8
rs2592555	11p13	PRR5L	C/T	102,760 (21)	0.27	0.93 (0.90–0.96)	8.7×10−7	116,551 (25)	6.78	AD	7
NOVEL LOCI
rs2038255	14q13.2	PPP2R3C	T/C	102,760 (21)	0.18	1.11 (1.07–1.14)	1.8×10−10	116,557 (25)	8.40
rs7127307	11q24.3	−/ETS1	C/T	103,066 (22)	0.47	0.93 (0.90–0.95)	3.9×10−10	116,855 (26)	9.08	SRA	14
rs7512552	1q21.2	C1orf51/MRPS21	T/C	102,762 (21)	0.49	0.93 (0.91–0.95)	9.1×10−10	116,544 (25)	6.91
rs6473227	8q21.13	MIR5708/ZBTB10	A/C	102,761 (21)	0.61	0.93 (0.91–0.95)	1.4×10−9	116,557 (25)	7.54	(AD),SRA,A+HF	9,14,53
rs6602364	10p15.1	IL15RA/IL2RA	G/C	103,065 (22)	0.45	1.08 (1.05–1.10)	1.5×10−9	116,855 (26)	7.86	(SRA)	14
rs10214237	5p13.2	IL7R/CAPSL	C/T	102,761 (21)	0.27	0.93 (0.90–0.95)	2.9×10−8	116,554 (25)	4.79
rs10199605	2p25.1	LINC00299/−	A/G	102,760 (21)	0.30	0.93 (0.90–0.95)	3.4×10−8	116,557 (25)	4.67	(SRA)	14
rs4643526	2p16.1	PUS10	A/G	103,066 (22)	0.19	1.09 (1.06–1.12)	3.5×10−8	107,425 (25)	6.31
rs12951971	17q21.2	STAT3	G/T	102,761 (21)	0.09	1.13 (1.08–1.17)	4.1×10−8	107,120 (24)	5.33
rs7625909	3p21.1	SFMBT1/RFT1	T/C	102,761 (21)	0.32	1.07 (1.05–1.10)	4.9×10−8	116,558 (25)	5.83
rs112111458	2p13.3	CD207/VAX2	G/A	102,760 (21)	0.13	0.91 (0.87–0.94)	1.4×10−7	116,553 (25)	6.57

in LD with known functional mutation rs2228145 (r2=0.86)

nearby SNP (rs6827756, bp position: 123184411) in LD (r2=0.97 in 1000genomes) showed similar association, log10BF=7.21, European fixed effects p-value 3×10−9

Nearest genes are the two flanking genes if intergenic (with the closer gene in bold, - indicates no gene within 250kb), single genes denote the variant is intronic.

at 1q21.2: variant is closest to LCE3A, but previously associated FLG is within 250kb, at 4q27: variant is within an intron of KIAA109, but previously associated IL2 is within 150kb, AD= atopic dermatitis, A=asthma, AS=allergic sensitization, SRA=self-reported allergy, AR=allergic rhinitis, A+HF=asthma and hayfever combined,

P-values and −log 10 Bayes Factors (BF) in bold indicate genome-wide significant results

EA/OA= effect allele/other allele, EAF = effect allele frequency, OR=odds ratio, CI=confidence interval, N= sample size, BF= bayes factor

Figure 1

Atopic dermatitis GWAS meta-analysis results

(A) Manhattan plot of European fixed effects meta-analysis. (B) Manhattan plot of the multi-ethnic MANTRA meta-analysis of all studies. Arrows mark variants not associated in the European-only analysis. (C) QQ plot of the European analysis - lambda=1.054.

Seven of the 21 established asthma loci[15-20], 7 of the 10 allergic sensitization loci[21], and 6 of 14 self-reported allergy loci[14] showed association with atopic dermatitis (p<0.05), all with consistent directions of effect, supporting common atopic mechanisms in atopic dermatitis and allergy (Supplementary Table 3). However, several studies used here contribute to multiple GWASs, which may bias this overlap. Nevertheless, a substantial proportion of the loci associated with other atopic conditions appear not to be strongly associated with atopic dermatitis.

Twenty-one of the 27 atopic dermatitis-associated loci have previously been implicated in other immune-mediated traits (Supplementary Table 4), most notably inflammatory bowel disease (IBD) and psoriasis. We therefore compared significant results from GWAS of IBD[22], psoriasis[23], ankylosing spondylitis[24], multiple sclerosis[25], rheumatoid arthritis[26] and type 1 diabetes[27]with results from our present study of atopic dermatitis. Of 163 established IBD risk variants, 39 reached p<0.05 for atopic dermatitis (Supplementary Table 5, 8.1 expected, p=2.4×10−16), 35 with the same direction of effect (sign test p<0.0001), consistent with the observational association between the two diseases[28-30]. Of the 36 known psoriasis variants, 15 reached p<0.05 for atopic dermatitis (Supplementary Table 6, 1.8 expected, p=6×10−11), 10 with the same direction of effect (sign test p=0.30). However, these conditions rarely clinically co-occur[31] and the most strongly associated genetic variants show opposite directions of effect[32]. Therefore our results, suggesting a more complex genetic relationship, might warrant further investigation. SNPs robustly associated with other auto-immune diseases were also more likely to be nominally associated with atopic dermatitis than expected by chance, but there was little evidence of any consistency in direction of effect (Supplementary Tables 7–10). These findings did not appear to be affected by contamination by common controls across studies. Analyses performed excluding common cohorts, yielded similar results (data not shown).

Conditional analysis showed evidence for secondary independent signals at 4 known atopic dermatitis loci (2q12.1, 4q27, 11p13, 5q31.1, Supplementary Table 11), one of which (5q31.1) has been previously reported[9]. In the epidermal differentiation complex (1q21.2–3, where FLG is located) the signals near MRPS21 (rs7512552) and IL6R (rs12730935 or the known functional mutation rs2228145) were independent from FLG, whereas the top signal near LCE3E (rs61813875) appears to be partially tagging the R501X FLG mutation (r2=0.49) and showed no significant residual association (P>0.05) after conditioning on the 4 most prevalent FLG mutations (Supplementary Tables 12&13).

To identify additional variants of biological relevance not reaching genome-wide significance, we applied gene-set enrichment analysis using Meta-Analysis Gene-set Enrichment of variaNT Associations (MAGENTA)[33] (Supplementary Table 14). A significant enrichment of 22 partially overlapping gene-sets (FDR<=0.01) related to innate immune signaling and T-cell polarization was observed (Supplementary Fig.5).

For replication, we selected the lead SNPs from the 11 novel loci, 9 candidate SNPs from the MAGENTA analysis (with p<10−5 mapping to gene-sets with FDR<0.05), and 3 SNPs representing potentially novel secondary signals. These were investigated in 18 studies (32,059 cases and 228,628 controls, Supplementary Table 1). Amongst the European studies, 11 of the 20 novel loci reached a Bonferroni-corrected threshold (α=0.0025) with 1-sided tests in a fixed effects analysis (Table 2). However, one of these showed evidence of heterogeneity (10p15.1, p=0.041) and was not significant in a random effects analysis (p=0.019, Supplementary Table 15). Two of the gene-set selected SNPs reached genome-wide significance in the combined analysis (2q37.1, 12q15). A random effects analysis of all replication cohorts (European and other ethnicities) show broadly consistent results (though only 6 reach genome-wide significance), with no clear population-specific effects (Supplementary Table 16 & Fig.6).

Table 2

Replication results for the novel genome-wide significant loci and loci identified in the MAGENTA gene-set enrichment analysis. Discovery, replication and combined results are shown.

					Discovery European			Replication European			Overall European - fixed effects			het	random effects p-values
Variant	Locus	Nearest Gene	EA/OA	EAF	N (studies)	OR (95% CI)	P-value	N (studies)	OR (95% CI)	P-value‡	N	OR (95% CI)	P-value	p-value	European	all studies
NOVEL GENOMEWIDE SIGNIFICANT LOCI
rs7512552	1q21.2	C1orf51/MRPS21	T/C	0.49	102762 (21)	0.93(0.91–0.95)	9.1×10−10	257019 (15)	0.98(0.96–0.99)	0.0048	359781	0.96(0.94–0.97)	5.41×10−9	0.002	1.5×10−7	1.3×10−7
rs10199605	2p25.1	LINC00299/−	A/G	0.30	102760 (21)	0.93(0.90–0.95)	3.4×10−8	256958 (15)	0.97(0.95–0.99)	0.0045	359718	0.96(0.94–0.97)	3.97×10−8	0.024	4.1×10−6	1.5×10−5
rs4643526	2p16.1	PUS10	A/G	0.19	103066 (22)	1.09(1.06–1.12)	3.5×10−8	257050 (14)	1.03(1.01–1.05)	0.0058	360116	1.05(1.03–1.07)	5.94×10−8	0.249	3.8×10−6	1.1×10−5
rs112111458	2p13.3	CD207/VAX2	G/A	0.13	102760 (21)	0.91(0.87–0.94)	1.4×10−7	257019 (15)	0.95(0.93–0.98)	7.95×10−4	359779	0.94(0.92–0.96)	9.38×10−9	0.076	4.4×10−6	1.6×10−7
rs11923593*	3p21.1	SFMBT1/RFT1	G/A	0.32	102761 (21)	1.07(1.04–1.10)	9.7×10−8	257002 (15)	1.01(0.99–1.03)	0.2600	359763	1.03(1.01–1.05)	1.30×10−4	0.081	8.7×10−5	1.7×10−5
rs10214237	5p13.2	IL7R/CAPSL	C/T	0.27	102761 (21)	0.93(0.90–0.95)	2.9×10−8	257010 (15)	0.94(0.93–0.96)	6.71×10−8	359771	0.94(0.92–0.95)	2.86×10−14	0.773	2.9×10−14	1.5×10−10
rs6473227	8q21.13	MIR5708/ZBTB10	A/C	0.61	102761 (21)	0.93(0.91–0.95)	1.4×10−9	257006 (15)	0.95(0.93–0.97)	4.53×10−9	359767	0.94(0.93–0.95)	2.22×10−16	0.622	2.2×10−16	5.3×10−18
rs6602364	10p15.1	IL15RA/IL2RA	G/C	0.45	103065 (22)	1.08(1.05–1.10)	1.5×10−9	256993 (15)	1.03(1.01–1.05)	3.91×10−4	360058	1.05(1.03–1.06)	1.33×10−10	0.041	3.6×10−6	1.6×10−6
rs7127307	11q24.3	−/ETS1	C/T	0.47	103066 (22)	0.93(0.90–0.95)	3.9×10−9	257034 (15)	0.94(0.93–0.96)	2.51×10−10	360100	0.94(0.92–0.95)	1.48×10−18	0.935	1.5×10−18	1.0×10−20
rs2143950*	14q13.2	PPP2R3C	T/C	0.17	102762 (21)	1.10(1.07–1.14)	6.8×10−10	249940 (12)	1.07(1.04–1.09)	9.92×10−8	352702	1.08(1.06–1.10)	1.78×10−15	0.092	4.8×10−7	8.6×10−10
rs17881320*	17q21.2	STAT3	T/G	0.08	96796 (19)	1.12(1.07–1.17)	2.0×10−6	249949 (12)	1.05(1.02–1.09)	1.47×10−3	346745	1.08(1.05–1.11)	1.41×10−7	0.405	6.2×10−7	2.6×10−6
MAGENTA GENE-SET ENRICHMENT ANALYSIS LOCI
rs1057258	2q37.1	INPP5D	T/C	0.18	101012 (21)	0.94(0.91–0.97)	6.57×10−5	257030 (15)	0.94(0.92–0.96)	3.79×10−7	358042	0.94(0.92–0.96)	1.72×10−10	0.811	1.7×10−10	4.1×10−13
rs6872156	5q35.1	DUSP1	A/G	0.24	103066 (22)	0.93(0.91–0.96)	2.35×10−6	257047 (15)	0.97(0.95–0.99)	0.0055	360113	0.96(0.94–0.97)	8.08×10−7	0.340	1.8×10−6	2.7×10−7
rs7016497	8q24.3	PTK2	T/C	0.21	103066 (22)	0.94(0.91–0.97)	1.09×10−4	257040 (15)	0.98(0.95–1.00)	0.0290	360106	0.96(0.95–0.98)	9.37×10−5	0.757	9.4×10−5	1.4×10−6
rs2905493	11q12.2	CD6/CD5	T/C	0.01	89617 (15)	0.78(0.68–0.89)	2.63×10−4	254992 (13)	1.01(0.94–1.08)	0.6040	344609	0.95(0.90–1.02)	0.1432	0.150	0.419	0.098
rs1799986	12q13.3	LRP1(STAT6)†	T/C	0.15	99165 (20)	0.91(0.88–0.94)	1.14×10−7	257022 (15)	0.98(0.96–1.01)	0.1140	356187	0.96(0.94–0.98)	2.90×10−5	0.182	1.1×10−4	1.6×10−3
rs2227483	12q15	IL22(& IFNG)†	A/T	0.44	102762 (21)	0.94(0.92–0.97)	2.27×10−6	253446 (14)	0.94(0.92–0.96)	3.55×10−11	356208	0.94(0.93–0.96)	6.66×10−16	0.664	6.7×10−16	1.2×10−17
rs7146581	14q32.32	TRAF3	T/C	0.24	102760 (21)	0.95(0.92–0.97)	1.44×10−4	256971 (15)	0.96(0.94–0.98)	5.67×10−5	359731	0.95(0.94–0.97)	6.17×10−8	0.219	1.2×10−4	4.1×10−6
rs11657987	17q25.3	PGS1(SOCS3)†	T/G	0.49	100695 (21)	1.06(1.04–1.09)	1.07×10−6	257019 (15)	1.03(1.01–1.05)	1.04×10−3	357714	1.04(1.03–1.06)	5.09×10−8	0.235	9.7×10−5	6.0×10−5
rs77714197	19q13.11	CEBPA	T/C	0.02	87690 (14)	1.35(1.19–1.54)	3.31×10−6	256447 (14)	0.98(0.89–1.08)	0.6540	344137	1.10(1.02–1.18)	0.0139	0.048	0.102	0.116

rs11923593 replaces rs7625909 (r2=0.98), rs2143950 replaces rs2038255 (r2=0.94), rs17881320 replaces rs12951971 (r2=0.75) in the replication analysis

rs1799986 is within LRP1, but was selected in the MAGENTA analysis due to its proximity to STAT6. rs2227483 is with IL22, but was selected due to its proximity to both IL22 and IFNG. rs11657987 is within PGS1, but was selected due to its proximity to SOCS3

Replication p-values for a 1-sided test

Replication p-values in bold were considered significant (p<0.0025), overall p-values in bold are genome-wide significant, heterogeneity p-values<0.05 are in bold

EA/OA= effect allele/other allele, EAF = effect allele frequency, OR=odds ratio, CI=confidence interval, N= sample size, het=heterogeneity

All 3 secondary signals showed significant association in the replication-phase conditional analysis (Supplementary Table 11).

As a preliminary step towards understanding the functional underpinnings of the atopic dermatitis genetic associations, we established a ‘credible set’ of SNPs (all with strong association) for each locus as described in the online methods[34]. We reviewed these SNPs’ functional annotations in ENCODE Consortium and Roadmap Epigenomics Consortium data, evaluated expression quantitative trait locus (eQTL) effects in MuTHER[35], reviewed evidence of differential expression, and surveyed relevant mouse mutants (see Supplementary Note 2 and Tables 17–21). Regions of DNase hypersensitivity from the ENCODE and Roadmap data[36,37] were strongly enriched for atopic dermatitis association compared to the rest of the genome (Supplementary Fig.7 & Table 22), particularly in immune cells (Th0, Th1, Th17 p<0.0001), this enrichment was observed well below the genome-wide significance threshold, indicating the presence of additional undetected risk variants. We observed multiple cis-eQTLs (Bonferroni-corrected p<7×10−4) in lymphoblastoid cell lines (LCLs) or skin (Supplementary Tables 17&19). The most significant were two variants from the credible set at 2p13.3, which were strong eQTLs for CD207/langerin in skin (rs4852714 p=1.23×10−10, rs6723629 p=1.67×10−10, LD with lead SNP r2=0.56, D’=0.96, and r2=0.53, D’=0.93, respectively, 99% posterior probability that atopic dermatitis and eQTL signals colocalize). rs4852714 is also in an open-chromatin region with histone marks indicative of promoter/enhancer activity in LCLs (Supplementary Tables 18,19 & Fig.8). CD207 encodes an intracellular pattern recognition receptor expressed in subpopulations of dendritic cells, in particular epidermal Langerhans cells (LCs) which play a vital role in the induction of tolerance and direction of adaptive immune responses[38]. CD207 binds to carbohydrates present e.g. on microorganisms and exerts anti-viral/anti-fungal defense mechanisms[39]. Of note, atopic dermatitis is characterized by an increased susceptibility towards skin infection with pathogens such as Staphylococcus aureus, herpes simplex virus, and Malassezia species[40], and differences in langerin function might contribute to this dysregulated cutaneous immunity.

There is longstanding evidence that skin barrier defects and inappropriate immune responses to environmental antigens[1] contributes to atopic dermatitis. However, evidence for autoimmune mechanisms, in particular in the context of progression to the chronic phase, has only recently emerged[41]. Interestingly, the majority of our novel susceptibility loci harbor candidate genes with functional annotations related to autoimmunity. At 14q13.2, the lead SNP (rs2038255) is intronic to PPP2R3C (a protein phosphatase component regulating B-cell maturation and survival), the dysregulation of which has been associated with murine autoimmunity[42] and the signal colocalizes with a strong KIAA0391 eQTL signal (Supplementary Table 19). The lead 5p13.2 variant (rs10214237) is located 4kb downstream of the gene encoding the alpha-chain of the IL7 receptor (IL7R), which is a key mediator in T-cell-driven autoimmunity and inflammation[43]. Of interest, the credible set contains an IL7R missense variant (rs6897932, p=1.6×10−7, r2=0.94 with lead SNP), which displays the same effect direction with multiple sclerosis[44,45]. The risk allele leads to an enhanced bioavailability of IL7[46], which in mice causes severe dermatitis with intense pruritus and high IgE levels, i.e. atopic dermatitis-like features[47]. Likewise, as part of the autosomal-dominant hyper-IgE syndrome, rare dominant negative mutations in the gene encoding STAT3 (in which our lead 17q21.2 variant is intronic) cause severe dermatitis and high serum IgE levels, as well as recurrent S.aureus skin infections, which may be driven by impaired Th17 cell differentiation and effector function[48,49]. STAT3 might thus represent an example for risk gene/pathway shared between a complex trait and a related Mendelian condition[50,51], harboring highly penetrant severe effect rare mutations and common milder effect variants. At 8q21.13, the closest candidate gene is ZBTB10 encoding a zinc finger protein, which is a putative repressor of the Sp1, Sp3 and Sp4 transcription factors[52]. Variants in moderate LD (r2>0.7) with the lead variant for atopic dermatitis were previously associated with self-reported allergy[14] and a combined asthma and hay fever phenotype[53]. However, although not excluding ZBTB10 as the causal gene, the credible SNP set comprises a 47kb interval on the other side of a recombination peak (60cM/Mb). The variant most likely to be regulatory amongst this set, deletion rs5892724 (r2=0.82 with lead SNP), is located in open chromatin in several cell types including CD4+ helper T-cells, and affects a STAT3 binding site[49,54]. At 11q24.3 the most plausible candidate gene is ETS1, which encodes a transcription factor with a range of immune functions including Th17 and B-cell differentiation and function; ETS1-deficient mice display autoimmune features[55]. ETS1 appears to be additionally involved in keratinocyte differentiation and formation of the cornified envelope[56]. Additional variants identified through the gene-set approach implicate genes with cytokine signaling functions (INPP5D, TRAF3, SOCS3 and a cytokine cluster on 12q15).

In conclusion, we have identified 10 new loci robustly associated with atopic dermatitis in Europeans (6 of which also reach genome-wide significance in random effects analysis across studies of all ethnicities), bringing the total number of susceptibility loci to 31 (24 in Europeans), with evidence of secondary signals at 4 of these. Altogether, in the subset of European studies with clinically defined cases, previously established and newly identified variants explain approximately 12.3% and 2.6% of the variance in liability, respectively (Supplementary Table 23). All novel susceptibility loci are related to (auto-)immune regulation, in particular innate signaling and T-cell activation and specification, and there appears to be a substantial genetic overlap with other inflammatory and autoimmune diseases. Whilst not detracting from the importance of maintaining the skin barrier in the prevention and treatment of atopic dermatitis, our findings lend support to new therapeutic approaches targeted at immune modulation[57].

ONLINE METHODS

GWAS meta-analysis

We carried out genome-wide association (GWA) analysis for atopic dermatitis case/control status in 26 individual studies (Supplementary Table 1), comprising a total of 21,399 cases and 95,464 controls. The majority of these studies included individuals of only European ancestry (22 studies, 18,900 cases, 84,166 controls). We also included one study of Japanese individuals (RIKEN, 1,472 cases. 7,966 controls), one study of African American individuals (SAPPHIRE, 422 cases and 844 controls), one study of Latin American individuals (GALA II, 300 cases, 1,592 controls) and one study with individuals of mixed non-European ancestry (Generation R, 305 cases, 896 controls).

Each cohort separately imputed their genetic data to 1000 Genomes Project Phase 1 (the majority to the March 2012 release, Supplementary Table 1) and carried out GWA analysis across all imputed variants. Before meta-analysis we restricted each study to only those variants with minor allele frequency (MAF)>1% and moderate imputation quality score (Rsq>0.3 for variants imputed in MACH and proper info>0.4 for IMPUTE). For some cohorts additional quality control filters were applied (full methods for each study are available in Supplementary Note 1).

Meta-analysis was conducted for Europeans only in GWAMA (using fixed effects) and for all ethnicities combined in MANTRA[60]. Rather than imposing a fixed or random effects model, MANTRA accounts for the heterogeneity of effects between ethnicities by allowing the studies to cluster according to allele frequency profile (and hence population genetic similarity). To prevent very small European studies (with less precise estimates of the allele frequencies) from having undue weight in our analysis we fixed the Europeans to cluster together by using the European fixed effects results in the MANTRA analysis. Variants with p<5×10−8 in the European analysis were considered to be associated, as were any additional variants with (log10) Bayes Factor (BF)>6.1 (equivalent to p<5×10−8)[61] in the MANTRA analysis. Each locus is represented in the results table by the variant with the strongest evidence for association. Heterogeneity was assessed using the I2 statistic and Cochrane’s Q test. Meta-analysis results were also stratified according to ethnicity, method of case diagnosis and age of onset to explore sources of heterogeneity.

For the Epidermal-differentiation complex region (where the FLG gene is located and which has previously shown complex association results), we repeated the association tests (across the region between 150.2–154.5Mb on chromosome 1) conditioning on the four most common FLG variants (R501X, 2282del4, R2447X, S3247X) in the individual studies where these were available (10 studies, 20,384 individuals, Supplementary Table 12). These were meta-analyzed to identify whether there were any remaining independent association signals in this region.

Identification of independent secondary signals at associated loci

To identify secondary independent signals at each of the other associated loci, we carried out conditional analysis of the European meta-analysis results using GCTA[62], with the ALSPAC 1000 Genomes imputation (restricted to variants with MAF>1% and imputation quality proper info score>0.8) serving as the reference. The regions tested were +/−250kb surrounding the top hit at each locus. Locus specific significance thresholds were estimated by first calculating the effective number of tests in each 500kb region using Nyholt’s procedure[63] and the 1000 Genomes reference data (European). For each locus we estimated the new threshold for locus-wide error rate of 5% by dividing alpha (0.05) by the corresponding number of effective tests in that region (α-values are shown in Supplementary Table 11). For 4q27 we defined the region as +/−500kb as a known hit was just less than 500kb from the top SNP in our analysis. We conditioned each region on the top hit from our meta-analysis. Any variant that surpassed the locus-specific threshold was considered an independent secondary signal.

MAGENTA gene-set enrichment analysis

We tested our meta-analysis results for enriched gene-sets using MAGENTA[33]. This method assigns SNPs to genes based on genomic distance (SNPs are assigned if within 110kb upstream or 40kb downstream of each gene), and generates gene-based summary p-values. Subsequently, genes are assigned to gene-sets (using curated repositories including GO-data, Biocarta, PANTHER, KEGG, etc.) and each gene-set is assigned a p-value by comparing gene summary p-values to a null model where SNPs are drawn randomly 10,000 times (normalizing for the number of SNPs genotyped in each gene) and controlling for false discovery rate (FDR) at α=0.05. ~10,000 gene-sets were tested. As MAGENTA requires a p-value as input and to take account of the differing effects between populations, we re-analyzed our meta-analysis of all studies using a random effects model, to serve as input to the MAGENTA analysis. All genes in the HLA region (chr6:29710331–33150000) were removed from the analysis. In order to identify additional variants of interest to take forward to replication we examined any pathway with FDR<0.05. From these gene-sets we took forward to replication any additional loci with a genetic variant p<10−5.

Cross-phenotype comparisons

The NHGRI GWAS catalog[64] was mined for traits with reported associations at each of our genome-wide significant loci. To further investigate the genetic overlap between atopic dermatitis and auto-immune diseases, we took the genome-wide significant loci from recent GWAS of IBD[22] and psoriasis[23] ankylosing spondylitis[24], multiple sclerosis[25], rheumatoid arthritis[26] and type 1 diabetes[27] and extracted the atopic dermatitis results for these variants from our European discovery GWAS, noting whether the variant was associated (p<0.05) with atopic dermatitis (testing enrichment of overlap using the 2-sided binomial exact test) and carried the same or opposite direction of effect between the two traits (tested using the sign test).

Replication

The top SNP from the 11 novel associated regions (log10BF>6.1 or p<5×10−8) were taken forward to replication, along with 9 suggestively associated SNPs (p<10−5) that were in genes highlighted in the MAGENTA analysis as good candidates. In addition, we included any variants with evidence for being secondary independent signals at associated loci. Replication consisted of 18 studies (32,059 cases and 228,628 controls) with genome-wide imputed data available or custom genotyping (Supplementary Table 1). Studies of European ethnicity were combined in fixed effects meta-analyses in GWAMA. We also carried out random effects meta-analysis for the European studies to assess association for variants which showed evidence of heterogeneity (p<0.05). Significant association in the replication phase was determined by 1-sided p-values meeting a Bonferroni-corrected threshold (α=0.05/20=0.0025). Random effects meta-analyses of replication studies of all ethnicities were also carried out and forest plots examined for evidence of population-specific effects. For the replication of the three secondary signals, the secondary SNPs were tested for association after conditioning on the top SNP in each of the European cohorts. These results were then combined in fixed effects meta-analyses.

Credible sets

In order to assemble a sensible list of variants at each locus for functional look-ups, we constructed credible sets[34] that represent those SNPs most likely to be causal based on statistical evidence from the MANTRA analysis (or from the European analysis for the three variants that appeared to be European-specific). The European-only GWAS was repeated in MANTRA to generate BFs required for the credible set analysis. Bayes factors were used to calculate posterior probabilities for all SNPs in each region (+/−1Mb), the minimum set of SNPs that had a cumulative posterior probability of at least 95% made up each credible set. These sets can be interpreted similar to confidence intervals, in that assuming the association signal at a locus can be attributed to a single causal variant (and that the true causal variant is included in the analysis and has been well-imputed), the 95% credible set contains that causal variant with 95% probability. Given that a 1000 Genomes imputation analysis may not include the true causal variant or that the associations may be driven by more than one causal variant, we do not expect these credible sets to necessarily contain the causal variants at the suggested 95% probability. Nevertheless, they demonstrate in addition to the ‘top SNP’, which neighbouring variants also show strong association with atopic dermatitis and we find them useful in assessing the size of the regions of interest and for defining a set of variants to follow-up. As the posterior probabilities for the MAGENTA identified credible sets are extremely large (due to the weaker signals at these loci), we instead carried out functional look-ups for all variants with r2>=0.8 of the top hit for these loci.

Functional look-ups

For variants identified as part of a credible set, we carried out look-ups in the following functional data resources; (i) regulomeDB and Haploreg were mined for evidence of coding or regulatory function (these resources collate annotations [e.g., coding variation, regulatory chromatin marks, DNase I hypersensitivity, protein binding and motif alteration] from the ENCODE Consortium, the NIH Roadmap Epigenomics Mapping Consortium, and the literature over a wide range of tissues); (ii) cis-eQTL for skin or LCLs were identified from the MuTHER consortium[35] (with variants considered eQTLs if association with any transcript within 1Mb was p<7×10−4, corresponding to a bonferroni correction for 36 loci and 2 tissues); (iii) differential expression reported for implicated genes between uninvolved skin from cases and skin from controls[65] and between lesional and non-lesional skin in atopic dermatitis patients in a study deposited in the Gene Expression Omnibus (Accession=GDS4444)[66]; and (iv) mouse mutants of implicated genes were examined in the MGI database.

Colocalization of atopic dermatitis GWAS signals and eQTLs in the MuTHER data was investigated using the R package coloc[67]. All SNPs within 100kb of the lead atopic dermatitis SNP were included in the analysis and we report the posterior probabilities that the two signals colocalize in Supplementary Table 19.

To identify and visualize cell types implicated in atopic dermatitis pathogenesis, the tendency of disease associated loci to fall in cell-type specific regulatory DNase Hypersensitive Sites (DHS) (a proxy for accessible and/or regulatory DNA) was calculated for the full range of p-values, essentially as done by Maurano et al.[68]. This enrichment was computed for 168 cell types and cell lines in the ENCODE Roadmap repository[36]. Duplicates and directly redundant cell types were removed before analysis. One-sided p-values for enrichment were calculated from an empirical null distribution of loci overlap for DHS-regions, generated by 10,000 random permutations of overlapping an identical number of random loci as found at p<=1×10−10 with all DHS-regions for all cell- and tissue types.

Variance in liability explained

We estimated the proportion of variance in atopic dermatitis liability explained by the established and novel variants in a subset of studies that had clinically diagnosed cases (GENEVA/KORAF4/POPGEN, NCRC-ADC, GENUFADex-SHIP1, GENUFAD-SHIP2, GENEVA(replication), CECCS) using the method of So et al. (2011)[69].

Data access

Genome-wide results are available on request to the corresponding author, on condition of signing any Data Transfer Agreements required according the IRB-approved protocols of contributing studies.