Meta-analysis of genome-wide association studies identifies three new risk loci for atopic dermatitis.

Atopic dermatitis (AD) is a commonly occurring chronic skin disease with high heritability. Apart from filaggrin (FLG), the genes influencing atopic dermatitis are largely unknown. We conducted a genome-wide association meta-analysis of 5,606 affected individuals and 20,565 controls from 16 population-based cohorts and then examined the ten most strongly associated new susceptibility loci in an additional 5,419 affected individuals and 19,833 controls from 14 studies. Three SNPs reached genome-wide significance in the discovery and replication cohorts combined, including rs479844 upstream of OVOL1 (odds ratio (OR) = 0.88, P = 1.1 × 10(-13)) and rs2164983 near ACTL9 (OR = 1.16, P = 7.1 × 10(-9)), both of which are near genes that have been implicated in epidermal proliferation and differentiation, as well as rs2897442 in KIF3A within the cytokine cluster at 5q31.1 (OR = 1.11, P = 3.8 × 10(-8)). We also replicated association with the FLG locus and with two recently identified association signals at 11q13.5 (rs7927894; P = 0.008) and 20q13.33 (rs6010620; P = 0.002). Our results underline the importance of both epidermal barrier function and immune dysregulation in atopic dermatitis pathogenesis.

Atopic dermatitis (AD), or eczema, is one of the most common chronic inflammatory skin diseases with prevalence rates of up to 20% in children and 3% in adults. It commonly starts during infancy and frequently precedes or co-occurs with food allergy, asthma and rhinitis[1] . AD shows a broad spectrum of clinical manifestations and is characterized by dry skin, intense pruritus, and a typical age-related distribution of inflammatory lesions with frequent bacterial and viral superinfections[1]. Profound alterations in skin barrier function and immunologic abnormalities are considered key components affecting the development and severity of AD, but the exact cellular and molecular mechanisms remain incompletely understood[1] .

There is substantial evidence in support of a strong genetic component in AD; however, knowledge on the genetic susceptibility to AD is rather limited[2,3]. So far, only null mutations in the epidermal structural protein filaggrin gene (FLG) have been established as major risk factors[4,5] .

The only genome-wide association study (GWAS) on AD in European populations so far identified a novel susceptibility locus on 11q13.5, downstream of C11orf30[6]. A recent second GWAS in a Chinese Han population identified two novel loci, one of which also showed evidence for association in a German sample (rs6010620, 20q13.33)[7]. In a collaborative effort to unravel additional risk genes for AD, we conducted a well powered two-staged genome-wide association meta-analysis in The EArly Genetics and Lifecourse Epidemiology (EAGLE) Consortium.

In the discovery analysis of 5,606 AD cases and 20,565 controls from 16 population-based cohorts of European descent (Supplementary Tables 1,2) there was little evidence for population stratification at study level (λGC<=1.08) or at the meta-analysis level (λGC=1.02), but an excess of association signals beyond those expected by chance (Supplementary Figs.1,2).

SNPs from two regions reached genome-wide significance in the discovery meta-analysis (Fig.1; Supplementary Table 3): rs7000782 (8q21.13, ZBTB10, OR=1.14, p=1.6×10−8) and rs9050 (1q21.3, TCHH, OR=1.33, p=1.9×10−8). Given the proximity of rs9050 to the well-established AD susceptibility gene FLG[4,5], we evaluated whether the observed association was due to linkage disequilibrium (LD) with FLG mutations. Despite low correlation between rs9050 and the two most prevalent FLG mutations (in ALSPAC (The Avon Longitudinal Study of Parents and Children): r2=0.257 for R501X, r2=0.001 for 2282del4) and high levels of recombination (peak of 20cM/Mb at ~150.4Mb in HapMap) between the TCHH and FLG regions, in a meta-analysis across eight studies conditional on the two FLG mutations, rs9050 was no longer associated with AD (OR=0.98, p=0.88) (Supplementary Fig.3) and was therefore not investigated further. rs9050 might tag a far-reaching haplotype on which the FLG null mutations occur, but we cannot exclude that there are additional AD risk variants in this complex region.

Figure 1

Manhattan plot for the discovery genome-wide association meta-analysis of atopic dermatitis

after excluding all SNPs MAF<1% and Rsqr<0.3 or proper_info<0.4. λ=1.017. SNPs with p<1×10−5 are shown in red.

The 11q13.5 locus previously reported to be associated in the only other European GWAS on AD to date[6] was confirmed in our meta-analysis (rs7927894 p=0.008, OR=1.07, 95%CI 1.02-1.12) (Supplementary Fig.4). So was the variant rs6010620 reported in a recent Chinese GWAS[7] (p=0.002, OR=1.09, 95%CI 1.03-1.15).

Of the 15 loci reported to be associated with asthma or total serum IgE levels in a recent GWAS[8], two showed suggestive evidence for association with AD (IL13:rs1295686, p=0.0008 and rs20541, p=0.0007; STAT6:rs167769 p=0.0379) (Supplementary Table 4).

After excluding the rs9050 SNP, we attempted to replicate the remaining 10 most strongly associated loci (P<10−5 in discovery, Table 1; Supplementary Table 3; Fig.2; Supplementary Fig.5) in 5,419 cases and 19,833 controls from 14 studies (Supplementary Tables 1,2). Three of the ten SNPs showed significant association after conservative Bonferroni correction (p<0.05/10=0.005) in the replication meta-analysis (and same direction of effect as the discovery meta-analysis): rs479844 near OVOL1, rs2164983 near ACTL9, and rs2897442 in intron 8 of KIF3A (Table 1; Fig. 2). All three SNPs reached genome-wide significance in the combined meta-analysis of discovery and replication sets: rs479844 (p=1.1×10−13, OR=0.88), rs2164983 (p=7.1×10−9, OR=1.16) and rs2897442 (p=3.8×10−8, OR=1.11). In contrast, rs7000782, which had reached genome-wide significance in the discovery analysis, showed no evidence of association in replication (p=0.296). There was no evidence of interactions between the three replicated SNPs (Supplementary Table 5).

Table 1

Discovery and replication results of the loci associated with atopic dermatitis

Results are for the fixed effect inverse-variance meta-analysis, with genomic control applied to the individual studies in the discovery meta-analysis. Stage I denotes discovery, II denotes replication and I+II denotes the combined analysis. The heterogeneity p-value (het pvalue), testing for overall heterogeneity between all discovery and replication studies was generated using Cochran’s Q-test for heterogeneity. All OR (odds ratios) are given with the minor allele representing the effect allele. CI denotes the confidence interval

Chr	SNP	Position(bp)	Gene	Effectallele	Otherallele	Effectallelefreq	Stage	N	OR (95% CI)	pvalue	hetpvalue
11	rs479844	65308533	OVOL1	A	G	0.44	I	26,151	0.89 (0.85,0.93)	7.8E-07
							II	25,098	0.87 (0.83,0.92)	2.4E-08
							I+II	51,249	0.88 (0.85,0.91)	1.1E-13	0.23
19	rs2164983*	8650381	ACTL9	A	C	0.15	I	17,403	1.22 (1.13,1.32)	1.8E-07
							II	22,996	1.11 (1.04,1.19)	0.002
							I+II	40,399	1.16 (1.10,1.22)	7.1E-09	0.004
5	rs2897442	132076926	KIF3A	C	T	0.29	I	26,164	1.12 (1.07,1.18)	7.8E-06
							II	25,064	1.09 (1.04,1.15)	0.001
							I+II	51,228	1.11 (1.07,1.15)	3.8E-08	0.52

rs2164983 was not included in the HapMap release 21 and so was missing for some discovery cohorts. This SNP showed evidence of heterogeneity (p=0.004). The random effects combined (I+II) result for this SNP was OR=1.14 (95%CI 1.05, 1.24) p=0.001.

Figure 2

Forest plots for the association of (a) rs479844, (b) rs2164983 and (c) rs2894772 with atopic dermatitis

All OR are reported with the minor allele (shown in brackets) as the effect allele. *MoBa imputation quality score was ‘info’ from PLINK.

GENR=Generation R. rs2164983 was not included in the HapMap release 21 and so was missing for some discovery cohorts.

Black points indicate the Odds Ratios (ORs) and the horizontal lines represent the 95% confidence intervals (CIs) for each study. Arrows are used to show where a CI extends beyond the range of the plot.

The sizes of the red and blue boxes indicate the relative weight of each study (using inverse variance weighting). Blue boxes indicate SNPs that were imputed and red boxes indicate SNPs on the genome-wide genotyping chip for the discovery cohorts and either on the genome-wide genotyping chip or individually genotyped for the replication cohorts. Only Health2006, KORA/GENEVA and NFBC86’ underwent individual SNP genotyping.

The subtotals (for discovery and replication) and overall ORs and CIs are indicated by the centre and range of the diamonds.

rs479844 (at 11q13.1) is located <3kb upstream of OVOL1. The pattern of LD is complex at this locus, but there is low recombination between rs479844 and this gene in Europeans (Supplementary Fig.2). OVOL1 belongs to a highly conserved family of genes involved in the regulation of the development and differentiation of epithelial tissues and germ cells[9-11] . It acts as a c-Myc repressor in keratinocytes, is activated by the β-catenin-LEF1 complex during epidermal differentiation, and represents a downstream target of Wg/Wnt and TGF-β/BMP7-Smad4 developmental signaling pathways[10,12,13]. Apart from their role in the organogenesis of skin and skin appendages[14,15], these pathways are also implicated in the postnatal regulation of epidermal proliferation and differentiation[16-18]. Disruption of OVOL1 in mice leads to keratinocyte hyperproliferation, hair shaft abnormalities, kidney cysts, and defective spermatogenesis[10,11]. In addition, OVOL1 regulates loricrin expression thereby preventing premature terminal differentiation[10]. Thus, it might be speculated whether variation at this locus influences epidermal proliferation and/or differentiation, which is known to be disturbed in AD. Analysis of transcript levels of all genes within 500 kb of rs479844 (OVOL1) in EBV-transformed lymphoblastoid cell lines (LCLs) from 949 ALSPAC individuals revealed a significant association between rs479844 and a nearby hypothetical protein DKFZp761E198 (p=7×10−5). Likewise, analysis of SNP-transcript pairs in the MuTHER (Multiple Tissue Human Expression Resource) skin genome-wide expression quantitative trait loci (eQTL) pilot database of 160 samples[19] provided suggestive evidence for an association in the same direction with DKFZp761E198 in one of the twin sets (Supplementary Fig.6). Further investigations are needed to clarify if the causal variant(s) at this locus exerts its effect through this transcript.

rs2164983 (at 19p13.2) is located in an intergenic region 70kb upstream of ADAMTS10 and 18kb downstream of ACTL9 (encoding a hypothetical protein). ADAMTS are a group of complex secreted zinc-dependent metalloproteinases, which bind to and cleave extracellular matrix components, and are involved in connective tissue remodelling and extracellular matrix turnover[20,21]. Actin proteins have well-characterized cytoskeletal functions, are important for the maintenance of epithelial morphology and cell migration, and have also been implicated in nuclear activities[22-24]. The low recombination between rs2164983 and ACTL9 and recombination peak between rs2164983 and ADAMTS10 in CEU HapMap individuals (Supplementary Fig.2) suggests the functional variant may be located within the ACTL9 region. There was no evidence for association between this SNP and any expression level of genes within 500kb in the ALSPAC LCL eQTL analysis, nor the MuTHER skin eQTL data (Supplementary Fig.6).

rs2897442 is located in intron 8 of KIF3A, which encodes a subunit of kinesin-II complex, required for the assembly of primary cilia, essential for Hedgehog signaling and implicated in β-catenin-dependent Wnt signaling to induce expression of a variety of genes that influence proliferation and apoptosis[25,26]. Of note, KIF3A is located in the 5q31 region, which is characterized by a complex LD pattern and contains a cluster of cytokine and immune-related genes, and has been linked to several autoimmune or inflammatory diseases, including psoriasis[27,28], Crohn’s disease[29,30], and asthma

[8,29,31] (Supplementary Table 4). In particular, distinct functional IL13 variants have been associated with asthma susceptibility[32]. Although rs2897442 is within the KIF3A gene, there is little recombination between this region and IL4 (interleukin 4). But there does appear to be a recombination peak between this region and IL13 (Supplementary Fig.7a). However, a secondary signal also appears to be present in the IL13/RAD50 region, and when conditioning on rs2897442 in our discovery meta-analysis, the signal in the IL13/RAD50 region remains, providing evidence of two independent signals (Supplementary Fig.7b). In an attempt to refine the association at this locus further, we analysed Immunochip data from 1,553 German AD cases and 3,640 population controls, 767 and 983 of which were part of the replication stage. The Immunochip is a custom content Illumina iSelect array focusing on autoimmune disorders, and offers an increased resolution at 5q31. In the population tested, the strongest signal was seen for the IL13 SNP rs848 (p=1.93×10−10), which is in high LD with the functional IL13 variant rs20541 (r2=0.979, D’=0.995). Further significant signals were observed for a cluster of tightly linked variants in IL4 (lead SNP rs66913936, p= 2.58×10−8) and KIF3A (rs2897442, p=8.84×10−7) (Supplementary Tables 6,7; Supplementary Fig.8). While rs2897442 showed only weak LD with rs848 (r2=0.160, D’=0.483), it was strongly correlated with rs66913936 (r2=0.858, D’=0.982). Likewise, pair-wise genotype-conditioned analyses showed that the significant association of rs2897442 with AD was abolished upon conditioning on rs66913936, whereas there was a remaining signal after conditioning on rs848 (Supplementary Tables 6,7). Analysis of LCL expression levels of all genes within 500kb of rs2897442 in ALSPAC revealed a modest association between rs2897442 and IL13 transcript levels (p=2.7×10−3). No associations with any transcript levels within 500kb of the proxy variant rs2299009 (r2=1) were seen in the MuTHER skin eQTL data (Supplementary Fig.6). However, this does not exclude a regulatory effect in another tissue or physiological state, involvement of causative variants in LD with these SNPs in long-range control of more distant genes[33], or different functional effects such as alternative splicing.

It is well known that genes that participate in the same pathway tend to be adjacent in the human genome and coordinately regulated[34]. Thus, our results and previous findings suggest that there are distinct effects at this locus, which might be part of a regulatory block. Further efforts including detailed sequencing and functional exploration are necessary to fully explore this locus.

Variants rs2164983, rs1327914 and rs10983837 showed evidence of heterogeneity in the meta-analysis (p<0.01). The overall random effects results for these variants were OR=1.14 (95%CI 1.05 −1.24), p=0.001; OR=1.06 (95%CI 1.00 - 1.13), p=0.058; and OR=1.11 (95%CI 0.98 - 1.20) p=0.155, respectively. Stratified analysis showed that the effects of rs2164983 and rs1327914 were stronger in the childhood AD cohorts (OR=1.23, p=2.9×10−9; OR=1.12, p=2.5×10−4) as compared to those studies that included AD cases of any age (OR=1.17, p=0.002; OR=1.02, p=0.584, p-value for the differences p=0.031 and p=0.028, respectively) (Supplementary Fig.9). This did not fully explain the heterogeneity for rs2164983 (in the childhood only cohorts the p-value for heterogeneity was still p<0.01). COPSAC (Copenhagen Studies on Asthma in Children) is noticeably in the opposite direction and excluding this study gives a heterogeneity p-value of 0.069 (OR=1.17, p=8.1×10−10). However, COPSAC is diagnostically and demographically comparable to the other cohorts and so there is no obvious reason why this cohort should give such a different result. Neither stratification by age of diagnosis nor whether a physician’s diagnosis was a case criterion explained the heterogeneity observed for rs10983837. Stratified analyses also indicated a stronger effect of rs2897442 in studies with a more stringent definition of AD (reported physician’s diagnosis) (OR=1.14, p=7.0×10−9) as compared to studies where AD was defined as self-reported history of AD only (OR=1.05, p=0.119) (Supplementary Fig.9). These observations underline the importance of careful phenotyping and support the claim of distinct disease entities rather than one illness as is reflected by current rather broad and inclusive concepts of AD. It is anticipated that the results of molecular studies will enable a more precise classification of AD.

In summary, in this large-scale GWAS on 11,025 AD cases and 40,398 controls we have identified and replicated two novel AD risk loci near genes which have annotations that suggest a role in epidermal proliferation and differentiation, supporting the importance of abnormalities in skin barrier function in the pathobiology of AD. In addition, we observed a genome-wide significant association signal from within the cytokine cluster on 5q31.1, this appeared to be due to two distinct signals, one centered on RAD50/IL13 and the other on IL4/KIF3A, both of which showed moderate association with IL13 expression. We further observed a signal in the epidermal differentiation complex, representing the FLG locus, and replicated the 11q13.5 variant identified in the only other (smaller) European GWAS on AD published to date. Our results are consistent with the hypothesis that AD is caused by both epidermal barrier abnormalities and immunological features. Further studies are needed to identify the causal variants at these loci and to understand the mechanisms through which they confer AD risk.