High-density genotyping of immune loci in Kawasaki disease and IVIG treatment response in European-American case-parent trio study.

Kawasaki disease (KD) is a diffuse and acute small-vessel vasculitis observed in children, and has genetic and autoimmune components. We genotyped 112 case-parent trios of European decent (confirmed by ancestry informative markers) using the immunoChip array, and performed association analyses with susceptibility to KD and intravenous immunoglobulin (IVIG) non-response. KD susceptibility was assessed using the transmission disequilibrium test, whereas IVIG non-response was evaluated using multivariable logistic regression analysis. We replicated single-nucleotide polymorphisms (SNPs) in three gene regions (FCGR, CD40/CDH22 and HLA-DQB2/HLA-DOB) that have been previously associated with KD and provide support to other findings of several novel SNPs in genes with a potential pathway in KD pathogenesis. SNP rs838143 in the 3'-untranslated region of the FUT1 gene (2.7 × 10(-5)) and rs9847915 in the intergenic region of LOC730109 | BRD7P2 (6.81 × 10(-7)) were the top hits for KD susceptibility in additive and dominant models, respectively. The top hits for IVIG responsiveness were rs1200332 in the intergenic region of BAZ1A | C14orf19 (1.4 × 10(-4)) and rs4889606 in the intron of the STX1B gene (6.95 × 10(-5)) in additive and dominant models, respectively. Our study suggests that genes and biological pathways involved in autoimmune diseases have an important role in the pathogenesis of KD and IVIG response mechanism.

Introduction

Kawasaki disease (KD) is the leading cause of acquired heart disease in children in the United States (U.S.). It is a potentially life-threatening acute vasculitis with diffuse involvement of multiple organ systems[1] and has a predilection for the coronary arteries. In the U.S. about 5500 cases were estimated in 2009 by passive surveillance alone;[2] and based on the system dynamics modeling simulations, there will be an average of 6200 new patients per year and a total of 161 776 individuals (half of them adults) with a history of KD by 2030.[3] The disease is primarily associated with the East Asian populations although there has been an increase in the incidence of KD among European-Americans over the last two decades.[4] The primary treatment of KD is the administration of intravenous immunoglobulins (IVIG) accompanied by aspirin.[5] Approximately 15 to 25% of the cases require more than a single dose of IVIG therapy to achieve complete resolution.[5] Furthermore, coronary artery aneurysms develop in 15-25% untreated cases compared to <5% among those treated, although a substantially higher percentage show coronary artery dilation.[5]

KD, also known as mucocutaneous lymph node syndrome (MCLS), is considered an autoimmune disease since walls of the blood vessels are inflamed throughout the body with damage to coronary arteries likely mediated by auto-antibodies to endothelial cells (ECs).[6] Parental KD is also considered a risk factor with the odds of KD among siblings and recurrent cases being 7 and 3 times higher among KD children and parents with history of KD, respectively compared to the general population,[7] suggesting a genetic etiology. However, genetic studies, including genome-wide association studies (GWAS) are often inconsistent and the associations are often race or ethnicity specific. Furthermore, only few candidate gene studies, including ours have assessed the gene-IVIG nonresponse association.[8-12] Our objective was to fine-map regions of immune-related genes using the Human ImmunoChip to validate the association of immune related genes with susceptibility to Kawasaki disease, search for novel SNPs, examine shared loci with other autoimmune diseases and perform genetic evaluation for IVIG treatment response among European American case-parent trios.

Results

The average age of the probands was 45.7 (± 32.3) months, with the average age for the responders being 44.8 (± 30.8) months and the non-responders being 46.6 (± 33.9) months. Overall, the study population consisted of 56.1% male probands, 48.1% males being responders and 64.2% being non-responders.

Susceptibility to Kawasaki disease

The Manhattan plots corresponding to the p-values for the SNPs in the transmission disequilibrium test (TDT) test are presented in Figure 1a for the additive model and Figure 1b for the dominant model. The Q-Q plots are illustrated in supplementary Figure S1a and S1b. The top 15 results from the FBAT analysis for the association between the SNPs and susceptibility to KD are presented in Table 1; upper panel for the additive model and lower panel for the dominant model (all results in supplementary Table 1 and Supplementary Figure S3 shows the LD plot for regions with multiple hits in chromosomes 1, 5 and 19). In the additive model, the top SNPs were predominantly located in chromosome 19, largely located in regions of galactoside 2-alpha-L-fucosyltransferase 1 (FUT1) gene and Ras interacting protein 1 (RASIP1) gene (Table1), but not in the region of inositol-trisphosphate 3-kinase C (ITPKC), another gene located on chromosome 19 that has been previously linked to KD. The topmost SNP (rs838143) and two other more common SNPs (rs838142 and rs4021), all in the untranslated region (UTR) of the FUT1 gene, were significantly over-transmitted from the parents to KD children (p=2.70×10-5, 4.80×10-5, 1.18×10-4, respectively).

Figure 1

Manhattan plot of the association between SNPs in the ImmunoChip and susceptibility to Kawasaki Disease using the transmission disequilibrium test (TDT) in a) additive model and b) dominant model

Table 1

Top 15 significant SNPs associated with susceptibility to Kawasaki disease using the a) additive model and b) dominant model in FBAT

CHR	SNP	Position	Gene Symbol	Gene Location	Risk* allele	MAF HAPMAP/1000 Genomes	p-Value
	Additive Model
19	rs838143	53943567	FUT1	3′ UTR	A	0.05	2.70 E-05
5	rs2054440	150365152	LOC134466 ∣ GPX3	Intergenic	G	0.42	3.00 E-05
19	rs838142	53943963	FUT1	3′ UTR	T	0.25	4.80 E-05
5	rs17111637	150361436	LOC134466 ∣ GPX3	Intergenic	A	0.35	7.30 E-05
3	rs12486347	70649298	LOC100128160 ∣ FOXP1	Intergenic	A	0.14	8.20 E-05
14	rs2273844	24173254	GZMB	5′ UTR	A	0.19	8.20 E-05
12	rs17539379	39414907	CNTN1	Intron	A	0.08	1.01 E-04
19	rs4021	53945073	FUT1	3′ UTR	G	0.25	1.18 E-04
19	rs463631	53926504	RASIP1	Intron	A	0.04	1.24 E-04
19	rs392775	53931669	RASIP1	Intron	T	0.04	1.24 E-04
5	rs6889028	55465663	ANKRD55 ∣ LOC727984	Intergenic	A	0.03	1.57 E-04
22	rs135058	42168027	MPPED1	Intron	T	0.05	1.62 E-04
5	rs2345000	150363456	LOC134466 ∣ GPX3	Intergenic	A	0.34	1.75 E-04
14	rs9919930	39686800	FBXO33 ∣ LRFN5	Intergenic	T	–	1.83 E-04
19	rs281384	53908738	MAMSTR	Intron	A	0.02	2.08 E-04
	Dominant Model
3	rs9847915	161279884	LOC730109 ∣ BRD7P2	Intergenic	A	0.29	6.81 E-07
11	rs118125355	35065834	PDHX ∣ CD44	Intergenic	–	0.01	6.82 E-07
1	rs4657067	159791329	FCGR3A ∣ FCGR2C	Intergenic	A	–	1.15 E-06
1	rs377982	159868496	FCGR3B ∣ FCGR2B	Intergenic	T	–	1.15 E-06
1	rs12735449	159785611	FCGR3A	Intron	T	–	1.36 E-06
1	rs76158727	159865216	FCGR3B	Intron	–	–	1.36 E-06
10	rs35980734	35301298	PARD3 ∣ CUL2	Intergenic	G	0.33	1.36 E-06
1	rs111974271	159847992	LOC1131982 ∣ FCGR3B	Intergenic	A	–	1.93 E-06
10	rs12251203	35479828	CREM	Intron	G	0.40	2.29 E-06
11	rs9651751	70956476	KRTAP5-10 ∣ KRTAP5-11	Intergenic	G	0.39	2.45 E-06
1	rs3964394	159803051	FCGR3A ∣ FCGR2C	Intergenic	–	–	2.73 E-06
1	rs143513265	151037507	LCE1D ∣ LCE1C	Intergenic	–	–	3.16 E-06
17	rs149259492	23182244	NOS2A ∣ LOC201229	Intergenic	–	0.35	3.16 E-06
16	rs762633	28515842	SULT1A2	5′ UTR	C	0.33	3.24 E-06
10	rs4325270	60185474	LOC728640 ∣ BICC1	Intergenic	A	0.55	4.30 E-06

Minor allele was always used as the risk allele in the analysis

In the dominant model, rs9847915 in the intergenic region between LOC730109 and Bromodomain containing 7 Pseudogene 2 (BRD7P2) on chromosome 3 was the topmost SNP significantly associated (p=6.81×10-7) with susceptibility to KD (Table 1, lower panel). However, association of SNPs was largely observed in genes on chromosome 1, in the region of FCGR gene family cluster.

IVIG response

The Manhattan plots that summarize the associations between each SNP and response to IVIG treatments are presented in Figure 2a for the additive model and Figure 2b for the dominant model (Q-Q plots in supplementary Figure S2a and S2b). The top 15 associated SNPs are presented in Table 2, upper panel for additive model and lower panel for dominant model (all results in supplementary Table 2 and Supplementary Figure S3 shows the LD plot for regions with multiple hits in chromosomes 6 and 16). In the additive model, intergenic variant rs1200332 between Bromodomain adjacent to zinc finger domain, 1A (BAZ1A) and C14orf19 genes on chromosome 14 was the topmost significantly associated SNP with response to IVIG treatment. The odds of not responding to IVIG treatment was 5 times higher in those with the minor allele T (OR-nonresponders: 5.49, p=1.40×10-4) compared to the wild-type C allele. On the other hand, individuals with a minor allele at SNPs rs6017164, rs72945401, rs7579420, rs7163190, rs3134926, rs3777914, rs1883137, and rs11641163 had higher odds of responding to IVIG and most of these responders had minor allele frequency (MAF) closer to the reported frequency for HapMap Caucasians (Table 2). Several SNPs on chromosome 6 were also found to be significantly associated with IVIG responsiveness in or around G-protein coupled receptor-6 (GPR6)/Wiskott-Aldrich syndrome protein (WASP)-family, member 1 (WASF1), protein reversionless 3-like (REV3L), neurogenic locus notch homolog protein 4 (NOTCH4)/ C6orf10, and tumor necrosis factor receptor family 3 interacting protein 2 (TRAF3IP2) genes.

Figure 2

Manhattan plot of the association between SNPs in the ImmunoChip and response to IVIG treatment using logistic regression in a) additive model and b) dominant model

Table 2

Adjusted Odds ratio (OR) for the top 15 SNPs associated with IVIG responsiveness using a) additive model and b) dominant model

CHR	SNP	Position	Gene Symbol	Gene Location	Minorallele	Odds Ratio-IVIG nonresponders	MAF-Responders	MAFNonresponders	MAFHAPMAP	p-Value
Additive Model
14	rs1200332	34436886	BAZ1A ∣ C14orf19	Intergenic	T	5.49	0.13	0.33	0.21	1.40 E-04
20	rs6017164	41815520	GTSF1L ∣ TOX2	Intergenic	A	0.25	0.45	0.32	0.45	1.95 E-04
6	rs4354185	110408354	GPR6	3′ UTR	A	4.46	0.25	0.43	0.48	2.74 E-04
13	rs7999399	88031506	SLITRK5 ∣ MIRH1	Intergenic	C	3.42	0.27	0.45	0.46	3.65 E-04
6	rs4317449	110408638	GPR6 ∣ WASF1	Intergenic	T	4.22	0.25	0.42	0.48	3.90 E-04
6	rs72945401	111833015	REV3L	Intron	G	0.02	0.14	0.02	0.12	5.82 E-04
5	rs10473594	22511857	CDH12	Intron	C	4.51	0.11	0.32	0.16	6.42 E-04
2	rs7579420	3938010	LOC728597 ∣ LOC727982	Intergenic	T	0.31	0.43	0.32	0.46	7.69 E-04
15	rs7163190	32812788	GOLGA8B ∣ GJD2	Intergenic	G	0.26	0.33	0.18	0.29	8.66 E-04
6	rs3134926	32308125	NOTCH4 ∣ C6orf10	Intergenic	C	0.27	0.40	0.25	0.20	9.12 E-04
6	rs3777914	112011860	TRAF3IP2	Intron	G	0.29	0.50	0.32	0.47	9.40 E-04
22	rs16995211	33781639	LARGE ∣ ISX	Intergenic	A	5.69	0.08	0.23	0.09	9.72 E-04
4	rs719379	40537658	APBB2	Intron	T	3.46	0.37	0.43	0.48	9.93 E-04
6	rs1883137	112013839	TRAF3IP2	Intron	G	0.29	0.47	0.34	0.46	9.98 E-04
16	rs11641163	79327374	CDYL2	Intron	A	0.24	0.43	0.37	0.47	1.07 E-03
Dominant Model
16	rs4889606	30918684	STX1B	Intron	G	0.13	0.49	0.27	0.35	6.95 E-05
16	rs9926533	30921680	STX1B	Intron	T	0.13	0.49	0.27	0.35	6.95 E-05
16	rs4889603	30889726	SETD1A	Intron	G	0.13	0.50	0.29	0.38	8.79 E-05
16	rs12917722	30855898	FBXL19	Intron	G	0.15	0.50	0.29	0.42	1.57 E-04
5	rs10473594	22511857	CDH12	Intron	C	5.82	0.11	0.32	0.16	1.87 E-04
14	rs1200332	34436886	BAZ1A ∣ C14orf19	Intergenic	T	5.56	0.13	0.33	0.21	2.25 E-04
16	rs8046707	30823734	CTF1 ∣ LOC283932	Intergenic	A	0.18	0.45	0.28	0.38	3.00 E-04
16	rs11649653	30825988	CTF1 ∣ LOC283932	Intergenic	G	0.18	0.45	0.28	0.38	3.00 E-04
6	rs4354185	110408354	GPR6	3′ UTR	A	6.22	0.25	0.43	0.48	3.09 E-04
6	rs3134926	32308125	NOTCH4 ∣ C6orf10	Intergenic	C	0.16	0.40	0.25	0.20	3.38 E-04
16	rs2288004	30961541	STX4 ∣ ZNF668	Intergenic	C	0.17	0.50	0.27	0.36	3.65 E-04
16	rs10871454	30955580	STX4	Intron	T	0.17	0.48	0.27	0.36	3.65 E-04
16	rs11150604	30944521	STX1B ∣ STX4	Intergenic	T	0.17	0.48	0.27	0.36	3.65 E-04
16	rs12448321	30955781	STX4	Intron	T	0.17	0.48	0.27	0.36	3.65 E-04
16	rs9939417	30960968	STX4 ∣ ZNF668	Intergenic	T	0.17	0.48	0.27	0.36	3.65 E-04

The top SNPs associated with response to IVIG treatment using the dominant model were predominantly located on chromosome 16 and several of them were in linkage disequilibrium (Table 2, lower panel). The top most significant SNP was rs4889606, located in the intronic region of the syntaxin 1B (STX1B) gene on chromosome 16 and in linkage disequilibrium with rs9926533. This SNP was associated with a lower odds of IVIG non-response among those carrying the minor allele G (OR-nonresponders: 0.13, p=6.95×10-5) compared to the wild-type “A” allele. Other SNPs located on chromosome 16 were found either in the intronic or intergenic regions of SET domain-containing protein 1A (SETD1A), F-box and Leucine-rich repeat protein 19 (FBXL19), cardiotrophin-1 (CTF1)/LOC283932, and syntaxin 1B and 4 (STX1B/STX4) or STX4/ zinc finger protein 668 (ZNF668) genes and all were associated with higher odds of response to IVIG treatment.

Discussions

In the present study, we evaluated the genetic association of susceptibility to KD and IVIG treatment response using the ImmunoChip in homogenous European-American case-parent trios. While the sample-size was relatively small, several SNPs from previous studies were replicated in our family-based study, strengthening our results. Furthermore, these replicated SNPs along with few of the top significant SNPs in our current study were found to be expression quantitative trait loci (eQTLs) for specific genes and are therefore more likely to play a functional role in the regulation of gene and expression and subsequent immune response to KD. Consequently, we discuss these SNPs in addition to some of the SNPs that were highly significantly in our study and have shown previous associations with autoimmune diseases.

Two SNPs, rs2857151 in the intergenic region of the human leucocyte antigen class II histocompatibility antigen, DQ beta 2 chain and DO beta (HLA-DQB2 and HLA-DOB), and rs4813003 downstream of TNF receptor superfamily member 5 (TNFRSF5) or B-cell associated molecule CD40 (CD40) gene that were previously reported in GWAS to be associated with KD were replicated (Table 3).[13] SNP rs2857151 was shown to be an eQTL for several of the HLA genes, namely, HLA-DOB, -DRB1, -DRB5, -DQA1, -DQA2, and transporter 2, ATP-binding cassette, sub-family B (TAP2).[14] SNP rs4813003 was shown to be an eQTL for the CD40 gene.[14, 15] SNP rs2857151 in combination with rs4813003 and another locus in a 3-locus model were previously shown to confer the highest risk for KD whereas the same SNP with rs1801274 was shown to confer the highest risk for CALs.[13] The compound effect was hypothesized to provide evidence for the additive effect of these loci on KD. SNPs within the TAP2 gene have also been associated with autoimmune diseases like rheumatoid arthritis (RA), ankylosing spondylitis and multiple sclerosis and the TAP2 proteins are known to play an important role in antigen presentation.[16-19] SNP rs4813003 is located in the intergenic region between the CD40 and the cadherin 22 (CDH22) genes. CD40 is a receptor protein which is part of the TNFR superfamily and is expressed on B cells, monocytes, dendritic cells, ECs and epithelial cells.[20] It interacts with the CD40 ligand on the T-helper cells and mediates immune and inflammatory responses. CD40 was also recently observed to facilitate damage to ECs through an inflammatory response mounted via activated platelets and the expression of adhesion molecules.[20] On the other hand, CDH22 is a type of cadherin, part of a superfamily of transmembrane proteins, and is known to mediate cell-to-cell adhesion in a calcium dependent manner. Altered expression of CDH22 has been previously associated with cell proliferation, invasiveness and metastasis in certain cancers.[21]

Table 3

Replication of previous genome-wide associations with susceptibility to Kawasaki Disease

CHR	SNP	Position	Gene Symbol	Gene Location	MAF HAPMAP	Model	p-Value
6	rs2857151	32871492	HLA-DQB2 ∣ HLA-DOB	Intergenic	0.44	Dominant	0.04
20	rs4813003	44196691	CD40 ∣ CDH22	Intergenic	0.16	Dominant	0.04

Previously, we had independently shown that rs1801274 in FCGR2A gene is associated with KD susceptibility in European-Americans,[9, 11] confirming similar findings from a GWAS study.[22] Multiple variants in the FCGR gene family cluster within the ImmunoChip were also associated with KD susceptibility in our TDT analysis with a dominant inheritance, with SNPs in both the intronic and intergenic regions associated with KD susceptibility. This gene region is consistent with our previous studies of KD.[9-11] These findings are encouraging, specifically since results from several GWAS studies have been replicated in only few studies using independent cohorts.[23, 24] Polymorphisms and copy number variation of FCGR genes are also associated with susceptibility for autoimmune diseases such as RA and systemic lupus erythematosus (SLE). Accordingly, an intersection could potentially exist between KD and other autoimmune disorders.

While the function of several genes in Table 1 is unknown or unconfirmed, few (FUT1, RASIP1, and BRD7P2) are noteworthy and could have potential biological mechanism in KD pathway. The SNPs significantly associated with KD in the additive model were located in the 3′UTR region of the FUT1 gene, the intronic region of the RASIP1 gene and the intergenic region of the BRD7P2 gene. FUT1 encodes the galactoside 2-alpha-L-fucosyltransferase 1 enzyme that is involved in the formation of the H antigen on erythroid cells, regulating the expression of ABO antigens.[25, 26] FUT1 plays a major role in the pathogenesis of RA and is involved in angiogenesis, cell adhesion and cell proliferation of various inflammatory cells.[27] This role is further supported by experiments showing inactivation of FUT1 using short-interfering RNA where reduction in cell proliferation and tumor growth were observed as well as migration of ECs and integration with the extracellular matrix (ECM).[28-30] Macrophages are also reported to modulate the expression of FUT1 and FUT2 through the interleukin-6 (IL-6) signaling pathway in epithelial cells of the endometrium to enable implantation of the embryo and increased levels of macrophages are potentially associated with inflammatory response.[31] The three SNPs from our study are located in the 3′ UTR of FUT1, and interestingly MAF for all three was higher in HapMap Asians compared to Caucasians but the same allele was the wild-type allele in Yoruban Africans. While the functions of these SNPs are not known, future studies could evaluate the potential binding sites of miRNAs in the 3′ UTR in relation to these variants. The RASIP1 gene encodes the RASIP1 protein that has been associated with tubulogenesis and angiogenesis in as early as the embryonic stage and is specific to the endothelium.[32] Lumen formation via RASIP1 occurs through the regulation of the Rho family of monomeric G proteins, requiring the suppression of RhoA and activation of Cdc42 and Rac1 proteins.[33] RASIP1 deficient vessels are shown to collapse as EC adhesion to the ECM is lacking.[33] RASIP1 variants have been associated with early microvascular disease, which is a predictor of incident cardiovascular disease, and Crohn's disease in conjunction with FUT2.[34, 35] Thus it may play a more important role in the pathogenesis of both KD and CALs. SNP rs9847915 in the BRD7P2 gene was associated with susceptibility to KD in our study; the MAF for this SNP is lowest among Africans but similar for Asians and Caucasians from HAPMAP. This gene has been associated with phenotypes such as Behcet's syndrome which also involves systemic vasculitis of small vessels with an auto-immune component.[36]

To date, IVIG responsiveness has not yet been evaluated using large-scale genotyping arrays. Even with a smaller sample-size, several genes (BAZ1A, GPR6, TRAF3IP2, REV3L, STX1B, STX4, and CTF1), listed in Table 2, are noteworthy as working mechanisms of IVIG therapy have yet to be fully elucidated. Additionally, SNP rs1200332 was shown to be an eQTL for chromosome 14 open reading frames 10 and 24 (C14orf10 and C14orf24) and SNP rs10871454 as an eQTL for ZNF668, STX4 and SETD1A among others.[14]

In the additive model, SNPs located in the intergenic region of the BAZ1A gene, the 3′UTR and intergenic regions of the GPR6 gene and the intronic regions of the TRAF3IP2 gene were significantly associated with IVIG responsiveness. The SNP rs1200332 in the intergenic region of BAZ1A and C14orf19 genes has been associated with Crohn's disease, an autoimmune disease, in Puerto Ricans.[37] The MAF of rs1200332 is higher among Africans but lower and similar across Caucasians and Asians in HAPMAP. The BAZ1A gene encodes a subunit of the adenosine triphosphate (ATP)-dependent chromatin assembly factor (ACF) which is part of the chromatic remodeling complex and is involved in DNA damage response and repair.[38] The GPR6 gene encodes the GPR6 protein that is expressed on vascular smooth muscles cells and ECs and is an orphan GPCR modulated by sphingosine-1 phosphate (S1P), an important regulator in the vascular system.[39] The TRAF3IP2 gene codes for the TRAF3-interacting protein 2, also known as activator 1 (Act1), that acts as a signaling adaptor in IL-17 mediated cellular immune responses through nuclear factor κ-B (NF- κB) activation or as a negative regulator by inhibiting CD40 and B-cell activating factor receptor (BAFFR) mediated signaling.[40] TRAF3IP2 has also been associated with auto-immune disorders like psoriasis and SLE and is speculated to bring about loss of the negative regulatory role of Act1 in SLE resulting in B-cell hyperplasia and enhanced survival of autoreactive B-cells.[40]

In the dominant model, most SNPs associated were in chromosome 16 in the intronic and intergenic regions of the STX gene family; the top SNP rs4889606 in the STX1B gene having the lowest MAF among Africans but the same minor allele being the wild-type allele among Asians from HAPMAP. In general, the interaction of syntaxin with synaptotagmin results in fusion of the vesicular and plasma membranes and subsequent exocytosis.[41] Smirnova et. al. speculated that syntaxin may be involved in the Lambert-Eaton syndrome, an autoimmune disorder associated with some small cell lung cancers resulting in antibodies against synaptotagmin and calcium channels.[42] Several SNPs within the intronic and intergenic region of the STX4 gene were associated with increased IVIG response. The STX4 gene encodes for STX4 which is also part of the syntaxin family of membrane integrated proteins.[43] STX4 plays an important role in the regulation of secretions from various immune cells including neutrophils, macrophages and eosinophils.[44] The addition of short interfering RNAs against STX4 has been shown to reduce IgE antibody exocytosis from plasma cells from several cell lines.[43] It was additionally shown that STX4 knock-down results in reduced secretion of tumor necrosis factor-α from macrophages.[44] SNPs in the intergenic region of the CTF1 gene were also significantly associated with IVIG responsiveness. The CTF1 gene encodes for the CTF1 protein which is part of the IL-6 family and is involved in cell proliferation, cell migration and collagen production. CTF1 is largely localized to the ECs and vascular smooth muscle cells with genetic variations in the gene associated with dilated cardiomyopathy or hypertension related left ventricular hypertrophy.[45-47] CTF1 has been shown to have both beneficial as well as detrimental effects based on which signaling pathway is induced since it has been associated with wound healing in myocardial infarction and prevention of apoptosis due to hypoxia as well as plaque instability and pathophysiology of atherosclerosis.[48, 49] The role of all these genes in IVIG response could potentially be mediated through various immunity signaling pathways.

There are several strengths to our study. Sample-size was relatively small in our study and thus the findings in our study are preliminary; however, considering that KD is rare (13.7 per 100,000 among whites < 5 years)[4] and only 15 to 25% do not respond to IVIG treatment, the difficulty in assembly of trios for this study should be appreciable. The trio-based family analysis using TDT is robust to factors like confounding and variance inflation which enabled us to confirm and replicate previous findings related to KD susceptibility despite the small sample size of our study. Additionally, we determined ethnicity based on ancestry informative markers (AIMs) and performed the analysis in confirmed European-Americans, accounting for population stratification and thus ensuring a homogenous population.

Our study is not without limitations. Of note, we included equal number of IVIG responders and non-responders for the TDT analysis; thus, some of the findings could be biased towards less frequent IVIG responders in the population, as they were preferentially included in the analyses. The ImmunoChip largely includes genes that were associated with several autoimmune disorders and was therefore ideal for testing genetic associations with KD but does not represent the whole spectrum of SNPs and genes that may be involved in immune response. As a result, the density of markers found on the ImmunoChip is not as high as opposed to fine mapping of gene regions and thus the SNPs associated in our study may be tagging one or more SNPs in LD and/or possibly the true causal variants. Lastly, we largely describe SNPs in the UTR and intronic regions although we acknowledge that variants from other parts of the region can be functional and related to distant genes.[50] And even though we include discussion of certain SNPs as eQTLs for some genes, we recognize that we were limited by the annotations provided by the ImmunoChip consortium and express caution in interpreting these results since certain eQTLs are reported in relation to specific cells and therefore may not be generalizable or applicable in the context of KD.

In the current study, we took a hypothesis-free agnostic approach to show that several SNPs in immune related genes included in the ImmunoChip are associated with susceptibility to KD and IVIG response. Some of these genetic associations have not yet been reported and can be potentially involved in KD pathogenesis while some of these genes are known to be in the pathway of other known autoimmune diseases. Future studies are therefore needed to replicate these findings and further fine-mapping is required to identify biologically causal SNPs as well as their functions in relation to either KD susceptibility or IVIG response mechanisms.

Methods

Study population

The current study included 112 case-parent trios of European-American descent from a retrospective cohort of KD. The study cohort has been described in detail previously.[11] Briefly, the probands (55 IVIG responders and 55 non-responders) were identified through clinical or hospital databases and then cross-referenced with the echocardiography databases of the respective heart centers. The probands were enrolled from children's hospitals at three locations: Seattle, WA; Oakland, CA; and Salt Lake City, UT. The study was approved by the Institutional Review Board at each of the participating institution and informed consent was obtained from the parents.

Clinical Outcomes

Kawasaki disease

Kawasaki disease was defined based on the guidelines recommended by the American Heart Association (AHA) and the American Academy of Pediatrics (AAP).[10, 51] The standard epidemiological criteria require the presence of fever for ≥ 5 days and at least 4 of 5 clinical features to make the diagnosis of KD. Participants in our study were included based on the initial criteria but we also included those with at least 2 clinical features accompanied by coronary artery involvement.

Response to IVIG

The treatment of Kawasaki disease included administration of an intravenous infusion of Immunoglobulin at 2 g/kg along with Aspirin within 10 days of the onset of fever. Nonresponse, or refractoriness to the initial IVIG therapy was characterized by the presence of persistent fever (≥ 38 °C) at ≥ 36 hours or recurrence of fever (≥ 38 °C) at ≥ 36 hours.[51]

DNA extraction, genotyping and quality control

As previously described, genomic DNA was extracted from blood or saliva samples obtained from the probands and their biological parents using the Versagene DNA purification kit (Gentra Systems, Minneapolis MN).[11] The PicoGreen assay for double stranded DNA was used for quantification and the final concentration was adjusted to 100ng/μl in Tris-EDTA.

Genotyping was performed using the ImmunoChip, an iSelect HD custom genotyping array (Illumina, San Diego CA), developed by the ImmunoChip Consortium.[52] The ImmunoChip contains 195,806 SNPs and 718 small insertions-deletions, of which 5,001 are non-synonymous coding, 1,926 are synonymous coding, 4,065 located in the untranslated region (UTR), and the rest fall either in the intronic, intergenic or undefined regions. The loci were selected on the basis of their associations with autoimmune or inflammatory diseases and ImmunoChip was designed using early 1000 Genome Pilot data (February 2010 release).[52, 53]

For genotype calling purposes, the data was analyzed using the Genome Studio Genotyping Module (Illumina, San Diego CA) using the NCBI build 36 (hg18) map. The normalized probe intensities were extracted for all samples that passed standard laboratory quality-control (QC) thresholds. Further QC included examining the SNPs for missingness by SNP and by individual, deviations from Hardy-Weinberg Equilibrium (HWE) at p-value <0.001 and minor allele frequencies < 0.05. Of a total 195,806 SNPs, there were 191,475 in the autosomal chromosomes. A total of 36 individuals were excluded as a result of missing > 10% genotypes data, of which 3 were probands and the rest were parents. Similarly, 16,220 SNPs were excluded since they were missing genotype data in > 10% of the sample population. Additionally, deviation from HWE was observed in 477 SNPs and 68,761 SNPs had minor allele frequency < 0.05 and therefore excluded from the analysis.

Statistical methods

TDT was employed to assess susceptibility to KD using family-based association test (FBAT), as described previously.[11] Briefly, parent-child trios were evaluated for differential transmission of alleles from heterozygous parents to the affected child. Single marker analysis was conducted to test the association under the null hypothesis of no linkage and no association using both additive and dominant genetic models. The direction and frequency of transmission is provided through the Z statistic and p-values.

The genetic association of response to IVIG treatment was assessed by implementing a case-control approach (responders versus non-responders) and performing logistic regression analyses adjusting for age and first two principal components, using both additive and dominant models in PLINK (v1.07).[54] European-American ancestry of all parents and probands was confirmed using 155 AIMs in EIGENSTRAT.[10, 11, 55] Quantile-quantile (Q-Q) and Manhattan plots, based on p-values, were created using SAS 9.3 (SAS Institute Inc., Cary, NC) to evaluate deviations of the observed test statistics distribution from the expected one and to visualize the association results, respectively.

Supplementary Figures

Quantile-quantile (QQ) plot of the association between SNPs and susceptibility to Kawasaki disease using a) additive model and b) dominant model

Quantile-quantile (QQ) plot of the association between SNPs and IVIG response to Kawasaki disease using a) additive model and b) dominant model

Locus-specific (-log10(p-value)) regional association plots for KD susceptibility (additive models in chromosomes 5 and 19 and dominant models in chromosome 1) and response to IVIG treatment (additive model in chrosmosome 6 and dominant model in chromosome 16) usin NCBI build 37, hg19. The diamond in purple indicates the most strongly associated signal. Estimated recombination rates are plotted in light blue graphical lines (labeled in right y-axis) to reflect the local LD structure from the HapMap population. The spectrum of colors indicates LD of each SNP with the most strongly associated signal (r2 values from HapMap where bright red indicates highly correlated, dark blue indicates weakly correlated and grey indicates missing r2 values).