Non-HLA type 1 diabetes genes modulate disease risk together with HLA-DQ and islet autoantibodies.

The possible interrelations between human leukocyte antigen (HLA)-DQ, non-HLA single-nucleotide polymorphisms (SNPs) and islet autoantibodies were investigated at clinical onset in 1-34-year-old type 1 diabetes (T1D) patients (n=305) and controls (n=203). Among the non-HLA SNPs reported by the Type 1 Diabetes Genetics Consortium, 24% were supported in this Swedish replication set including that the increased risk of minor PTPN22 allele and high-risk HLA was modified by GAD65 autoantibodies. The association between T1D and the minor AA+AC genotype in ERBB3 gene was stronger among IA-2 autoantibody-positive patients (comparison P=0.047). The association between T1D and the common insulin (AA) genotype was stronger among insulin autoantibody (IAA)-positive patients (comparison P=0.008). In contrast, the association between T1D and unidentified 26471 gene was stronger among IAA-negative (comparison P=0.049) and IA-2 autoantibody-negative (comparison P=0.052) patients. Finally, the association between IL2RA and T1D was stronger among IAA-positive than among IAA-negative patients (comparison P=0.028). These results suggest that the increased risk of T1D by non-HLA genes is often modified by both islet autoantibodies and HLA-DQ. The interactions between non-HLA genes, islet autoantibodies and HLA-DQ should be taken into account in T1D prediction studies as well as in prevention trials aimed at inducing immunological tolerance to islet autoantigens.

Introduction

Autoimmune (type 1) diabetes arises after a chronic prodrome of islet autoimmunity associated with a progressive loss of the pancreatic islet beta cells (for reviews see[1, 2]). The islet autoimmunity may be on-going months to years prior to the clinical diagnosis of diabetes [3-5]. The onset of hyperglycemia is the likely final step in a prolonged disease process that eventually results in insufficient beta-cell function to maintain normal blood glucose. The beta-cell autoimmunity that develops long before the clinical onset is becoming better understood through longitudinal investigations from birth of children at increased genetic risk for type 1 diabetes (T1D) [3, 5-8]. The apparent two-step pathogenesis of autoimmune (type 1) diabetes complicates the dissection of both the genetic etiology and the actual pathogenesis. The first step is perceived as an event(s) that triggers autoantibodies against one or several of the islet autoantigens, insulin, GAD65, IA-2 or ZnT8 [8-10]. The first autoantibody appears to be against either insulin, GAD65, or both [11]. Insulin autoantibodies tend to appear in HLA-DQ8 children with a peak during the first two years of life. GAD65 autoantibodies as the first autoantibody tend to appear in HLA-DQ2 children with a peak around 5 years of age [11]. The second step represents the time period from seroconversion to the clinical onset of diabetes. The length of the second step may be dependent not only on the number of islet autoantibodies [8, 12, 13] but also on HLA and non-HLA genes [14-16] as well as environmental factors [17].

HLA-DQ on the short arm of chromosome 6 is the most important gene for T1D [18, 19]. While certain DQ haplotypes are strongly associated with T1D, it remains to be fully established whether the principal association is with the risk for islet autoimmunity (step one), progress to clinical onset (step two), or both. It is important to distinguish between the association of HLADQ with the ability of an individual to mount an immune response to certain islet autoantigens and that between HLA-DQ and the loss of beta cells that would be consistent with the clinical diagnosis of diabetes [20, 21].

At the time of clinical diagnosis a major proportion of beta cells have been lost in a process involving macrophages or dendritic cells as potential antigen presenting cells as well as both T and B lymphocytes [1]. In the absence of reliable and reproducible cellular tests [22], the autoimmune pathogenesis and progression to clinical onset is best reflected by autoantibodies against insulin, GAD65, IA-2 and ZnT8 [23-26]. Also at the time of clinical diagnosis, autoantibodies against GAD65 are primarily associated with HLA-DQ2 [21], IA-2 with DR4 subtypes rather than with DQ8 [24, 27], insulin with DQ8 [21] and ZnT8 with DQ6.4 [9, 28]. We also observed that IA-2 autoantibodies are negatively associated with DQ2 [21, 28]. While the associations between HLADQ and islet autoantibodies may be mechanistically related to peptide presentations by DQ heterodimers (for reviews see [29, 30]), non-HLA genes also contribute to the risk for islet autoimmunity, progression to clinical onset of diabetes, or both [14, 31-33]. In a previous investigation of a Swedish population-based case-control study of 1-34 year old T1D patients and controls it was observed that both INS and DQ8 contributed to age-dependent risk for T1D [21]. The previously reported association between PTPN22 and T1D [34] was not only supported in our Swedish case-control study but also extended to demonstrate that the largest fold increase in risk of GAD65 autoantibody positive T1D associated with high risk PTPN22 allele was observed among those in the low risk HLA-DQ group [14].

During the course of investigating possible interrelations between HLA-DQ, non-HLA genes and islet autoantibodies [14, 21], the T1D Genetics Consortium (T1DGC) [35] completed Genome Wide Association Studies (GWAS) and reported at least 50 loci that conferred risk for T1D [15, 16]. An opportunity emerged from an invitation of the T1DGC to type as a replication set the Swedish cohort of patients and controls [21]. As all of the T1D patients and controls in this population-based case-control study stemmed from 1985-1989, thereby preceding the T1DGC effort, the first aim of the present study was to examine to what extent the T1D-associated non-HLA genes and loci reported by the TIDGC in much larger, and heterogeneous, datasets could be replicated in a smaller dataset from Sweden. The second aim was to find supporting evidence for our previous result that the increased risk of T1D associated with the minor (T) PTPN22 allele was modified by both HLA-DQ and autoantibodies against GAD65 [14]. The third aim was to estimate the associations between T1D and each of the non-HLA genes, stratified by HLA-DQ, as well as that between the non-HLA genes and autoantibodies against any of the islet autoantigens, GAD65, IA-2, or insulin as well as ICA, all measured within days after the clinical onset of diabetes [21].

Results

Type 1 diabetes and HLA as well as non-HLA genes

The T1DGC reported a total of 30 SNPs where the minor allele was related to increased risk for T1D. In the present data set we found supporting evidence for seven (DQA1, DQA2, PTPN22, C1QTNF6, HORMAD2, 730217 and ERBB3) genes of the minor allele non-HLA genes reported by the T1DGC (7/30, 23%; Table 1A). Among the total of 21 genes where the T1DGC reported the major allele to be associated with T1D, five (26471, INS, BACH2, IL2RA and SCAP2) were nominally statistically significantly associated with T1D in our study (5/21, 24%; Table 1B). In the analysis adjusted for HLA based on DQ genotypes, DQA1, PTPN22, C1QTNF6 were nominally statistically significantly associated with T1D in our study among the genes with the minor allele being the risk carrying allele for T1D (Supplementary Table 1A, Supplementary Figure 1A). Lastly, among the genes with the major allele associated with T1D, 26471, INS, BACH2, IL2RA were found to be nominally statistically significantly associated with T1D in our study (Supplementary Table 1B, Supplementary Figure 1B). The estimates based on our dataset for all 51 SNPs reported by the T1DGC are summarized in Tables 1A, 1B, Supplementary Tables 1A and 1B and Supplementary Figures 1A and 1B.

Table 1A

Summary of estimated OR (95% CI) of being diagnosed with type 1 diabetes for subjects with the risk carrying minor allele (including those with a heterozygous genotype mm+Mm) compared to subjects with the major genotype (MM), adjusting for age, sex, region of birth and HLA risk category (original). The number of patients and controls in the “risk carrying minor allele” columns count subjects with either the mm or Mm genotype. The number of patients and controls in the “major allele” column count subjects with the MM genotype. The SNPs are sorted by the p-value associated with the OR (95% CI). The corresponding results where the major allele carries risk are summarized in Table 1B. Lastly, the results adjusted for HLA risk categories based on DQ genotypes can be found in the Supplementary Tables 1A and 1B.

Original HLA risk categories
Gene	SNP	Risk carrying minor allele			Major allele			Total n for current SNP	OR (95% CI)	p-value
		Allele	Patients	Controls	Allele	Patients	Controls
HLA-DQA1	rs2187668	A	142	50	G	136	140	468	2.51 (1.55, 4.05)	1.72e-4
PTPN22	rs2476601	T	103	44	C	167	141	455	2.02 (1.24, 3.29)	0.005
C1QTNF6	rs229541	T	206	119	C	70	73	468	1.87 (1.17, 3.01)	0.009
HLA-DQA2	rs7454108	C	187	52	T	91	135	465	2.39 (1.18, 4.87)	0.016
HORMAD2	rs5753037	T	174	107	C	106	83	470	1.58 (1.01, 2.46)	0.043
730217	rs4900384	G	166	101	A	118	89	474	1.57 (1.01, 2.43)	0.046
ERBB3	rs2292239	A	178	109	C	101	81	469	1.58 (1.01, 2.48)	0.047
SH2B3	rs3184504	T	217	136	C	61	56	470	1.51 (0.92, 2.49)	0.106
442263	rs10499194	T	115	73	C	159	118	465	1.32 (0.85, 2.06)	0.220
CTRB	rs7202877	G	72	36	T	210	155	473	1.32 (0.78, 2.24)	0.294
TAGAP	rs1738074	G	228	151	A	44	41	464	1.27 (0.73, 2.22)	0.398
645481	rs10517086	A	137	86	G	144	102	469	1.20 (0.78, 1.84)	0.412
PTPN2	rs478582	C	205	138	T	74	52	469	1.22 (0.75, 1.98)	0.425
IL18RAP	rs917997	A	112	73	G	162	114	461	1.19 (0.76, 1.86)	0.442
IFIH1	rs1990760	A	238	157	G	32	26	453	1.29 (0.67, 2.50)	0.449
CCR5	rs333	DEL	13	11	INS	48	40	112	1.57 (0.49, 5.06)	0.451
SMARCE1	rs7221109	T	155	102	C	127	89	473	1.17 (0.76, 1.80)	0.479
UBASH3A	rs3788013	A	182	116	C	101	75	474	1.17 (0.75, 1.82)	0.494
IL2	rs2069762	G	122	78	T	160	113	473	1.15 (0.75, 1.77)	0.521
COBL	rs4948088	A	25	16	C	254	173	468	1.28 (0.60, 2.75)	0.523
C6orf173	rs9388489	G	194	129	A	85	59	467	1.14 (0.72, 1.82)	0.570
GAB3	rs2664170	G	114	81	A	165	107	467	1.11 (0.71, 1.73)	0.651
IL2RA	rs2104286	G	125	85	A	149	106	465	1.07 (0.69, 1.66)	0.758
IL21	rs6822844	T	91	60	G	188	130	469	1.06 (0.67, 1.68)	0.799
SIRPG	rs2281808	T	161	108	C	120	83	472	1.06 (0.68, 1.63)	0.805
C10orf59	rs10509540	C	135	88	T	148	105	476	1.05 (0.68, 1.62)	0.835
KIF5A	rs1678542	C	188	130	G	92	61	471	1.04 (0.66, 1.65)	0.860
IL10	rs3024505	T	80	59	C	213	137	489	1.04 (0.65, 1.66)	0.873
GSDML	rs2290400	A	200	130	G	82	57	469	1.02 (0.64, 1.64)	0.928
GLIS3	rs7020673	C	208	137	G	72	47	464	1.02 (0.62, 1.67)	0.950

Table 1B

Summary of estimated OR (95% CI) of being diagnosed with type 1 diabetes for subjects with the risk carrying major allele (MM) compared to subjects with the minor genotype (including those with a heterozygous genotype mm+Mm), adjusting for age, sex, region of birth and HLA risk category (original). The number of patients and controls in the “risk carrying major allele” columns count subjects with the MM genotype. The number of patients and controls in the “minor allele” column count subjects with either the mm or Mm genotype. The SNPs are sorted by the p-value associated with the OR (95% CI). The corresponding results where the minor allele carries risk are summarized in Table 1A. Lastly, the results adjusted for HLA risk categories based on DQ genotypes can be found in the Supplementary Tables 1A and 1B.

Original HLA risk category
Gene	SNP	Risk carrying major allele			Minor allele			Total n for current SNP	OR (95% CI)	p-value
		Allele	Patients	Controls	Allele	Patients	Controls
26471	rs4788084	G	108	42	A	171	147	468	2.33 (1.44, 3.77)	0.001
INS	rs689	A	195	96	T	82	93	466	2.00 (1.28, 3.11)	0.002
IL2RA	rs12722495	A	222	149	G	41	39	451	1.90 (1.09, 3.32)	0.025
SCAP2	rs7804356	T	195	116	C	90	74	475	1.65 (1.04, 2.61)	0.034
BACH2	rs11755527	C	56	27	G	218	160	461	1.89 (1.05, 3.43)	0.035
FHOD3	rs2644261	T	147	89	C	131	92	459	1.52 (0.97, 2.37)	0.067
RGS1	rs2816316	T	190	122	G	89	68	469	1.48 (0.93, 2.34)	0.096
DKFZp667F0711	rs947474	A	205	130	G	74	63	472	1.43 (0.90, 2.28)	0.127
CD69	rs4763879	G	76	41	A	136	92	345	1.34 (0.79, 2.29)	0.280
IL2	rs2069763	G	93	53	T	179	138	463	1.29 (0.81, 2.05)	0.292
CTLA4	rs3087243	C	120	74	T	151	116	461	1.26 (0.80, 1.98)	0.311
TMSB4X	rs5979785	T	186	112	C	90	72	460	1.26 (0.79, 2.00)	0.337
PRKD2	rs425105	A	212	142	G	67	46	467	1.20 (0.72, 2.01)	0.477
PTPN2	rs45450798	G	180	125	C	98	69	472	1.16 (0.74, 1.81)	0.527
CD226	rs763361	C	87	56	T	193	136	472	1.15 (0.72, 1.85)	0.550
ZFP36L1	rs1465788	G	143	93	A	134	96	466	1.08 (0.70, 1.66)	0.718
RBM17	rs11594656	T	185	121	A	94	69	469	1.07 (0.68, 1.69)	0.761
IL7R	rs6897932	C	166	110	T	113	79	468	1.07 (0.69, 1.67)	0.765
442263	rs6920220	G	155	106	A	122	86	469	1.07 (0.69, 1.64)	0.773
CLEC16A	rs12708716	A	138	97	G	125	89	449	1.06 (0.68, 1.64)	0.811
CTSH	rs3825932	C	116	84	T	162	109	471	1.01 (0.65, 1.57)	0.961

PTPN22, low risk HLA-DQ and GADA

We found that the minor (T) allele of PTPN22 contributed to an increase in risk of T1D, primarily in patients with neutral risk HLA-DQ (bottom panel of Table 2). This result, based on n=508 subjects, is consistent with our previous result for PTPN22 based on n=1240 subjects (top panel of Table 2) [14]. Both, our current dataset (n=508) and the dataset used for the previous PTPN22 analysis (n=1240) are subsets of the original matched case-control study (n=1708). This comparison between the prior [14] and current results for PTPN22 demonstrates that the reduced sample size available for typing of the remaining T1DGC identified SNPs was sufficient to find supporting evidence of a nominally statistically significant increase in T1D risk among carriers of the minor PTPN22 allele.

Table 2

Estimated odds ratios (OR) and 95% confidence intervals (95% CI) of being diagnosed with T1D for subjects with PTPN22(CT+TT) genotype compared with subjects with PTPN22(CC) genotype overall and stratified by HLA risk category adjusting for age, sex, and region of birth.

PTPN22 dataset (n=1240)[14]
HLA risk category(a)	PTPN22 WT (CC)		PTPN22 MUT (CT + TT)		Fold increase (MUT/WT)
	nP/nC	OR (95% CI)	nP/nC	OR (95% CI)	OR (95% CI)
All patients	444/404	1.00 (reference)	268/124	2.11 (1.57, 2.83)	2.11 (1.57, 2.83)
Low	44/194	1.00 (reference)	32/62	2.27 (1.32, 3.92)	2.27 (1.32, 3.92)
Neutral	110/143	3.52 (2.32, 5.34)	77/37	9.64 (5.74, 16.17)	2.74 (1.71, 4.38)
High	290/67	20.65 (13.44, 31.72)	159/25	29.26 (17.06, 50.19)	1.42 (0.86, 2.34)

(a) Low = {(DQB1*0602 or DQB1*0603), (DQB1*0301 but not DQB1*02, DQB1*0302, DQB1*0604)}, Neutral = {(DQA1*0501-B1*0201 but not DQA1*0301-B1*0302), Other genotypes (not appearing in low or high categories)}, High = {(DQA1*0501-B1*0201 and DQA1*0301-B1*0302), (DQA1*0301-B1*0302 but not DQA1*0501-B1*0201)}.

(b) n = 455 of these have rs2476601 (PTPN22 SNP) information ascertained by the T1DGC.

In addition, in the original data set, among GADA-positive patients, the OR was 2.83 (95%CI=2.00-3.99), whereas in GADA-negative patients, the OR was 1.41 (95%CI=0.98-2.04) (p for comparison=0.007). Furthermore, the OR of association between PTPN22 (CT+TT) and GADA-positive patients declined with increasing HLA-risk category from 6.12 to 1.54 (p=0.003); no such change was detected in GADA-negative patients (p=0.722) (p for comparison=0.001)[14]. In the current study, based on the analysis adjusted for age, sex, region of birth and HLA risk category the OR of association between PTPN22 (CT+TT) and GADA-positive T1D patients (OR=2.98 (1.69, 5.26); p=1.68e-4) was higher compared to that for GADA-negative T1D patients (OR=1.17 (0.62, 2.22)) (p=0.635) and this difference in association was nominally statistically significant (p for comparison=0.028) (Table 3, Figure 1).

Table 3

(Top) Estimated odds ratios (OR) and 95% confidence intervals (95% CI) of being diagnosed with T1D for subjects with the variant (minor) genotype compared to subjects with the major genotype and with low risk HLA haplotype (reference group), adjusting for age, sex and region of birth. The fold increase column summarizes the estimated OR (95% CI) of being diagnosed with T1D for carriers of the variant allele compared to non-carriers, for subjects in a given HLA risk category. The p-value = 0.005 for the model without interaction (top row), p-value = 0.016 for the model with interaction of PTPN22 and HLA risk category, and the p-value for trend = 0.705. (Bottom) Estimated OR and 95% CI of being diagnosed with autoantibody positive T1D (or autoantibody negative T1D) for subjects with the minor genotype compared to those with the major genotype, adjusted for age, sex, region of birth and HLA risk category.

HLA risk Category	PTPN22 (CC)		PTPN22 (CT+TT)		Fold increase OR (95% CI)
	nP/nC	OR (95% CI)	nP/nC	OR (95% CI)
	167/141	(reference)	103/44	2.02 (1.24, 3.29)	2.02 (1.24, 3.29)
Low	15/67	(reference)	8/21	1.61 (0.59, 4.41)	1.61 (0.59, 4.41)
Neutral	43/48	4.34 (2.13, 8.85)	37/13	13.3 (5.63, 31.43)	3.06 (1.43, 6.58)
High	109/26	19.76 (9.56, 40.83)	58/10	27.88 (11.41, 68.12)	1.41 (0.63, 3.16)

Low HLA risk = {(DQB1*0602 or DQB1*0603), (DQB1*0301 but not DQB1*02, DQB1*0302, DQB1*0604)}

Neutral HLA risk = {(DQA1*0501-B1*0201 but not DQA1*0301-B1*0302), Other genotypes (not appearing in low or high categories)}

High HLA risk = {(DQA1*0501-B1*0201 and DQA1*0301-B1*0302), (DQA1*0301-B1*0302 but not DQA1*0501-B1*0201)}

Figure 1

Summary of the estimated OR (dot) and 95% confidence intervals (whiskers) of the association between a given SNP and autoantibody-positive (-negative) T1D, adjusting for age, sex, region of birth and HLA risk category (original or one based on DQ genotypes). The red pairs indicate that the association between a given SNP and autoantibody-positive T1D was nominally statistically significantly different from the association between that SNP and autoantibody-negative T1D at 0.05 level (comparison p-value<0.05), the blue pairs indicate the same at 0.10 level (comparison p-value<0.10). The halos indicate that the OR is nominally statistically significantly different from one. Select SNPs are presented here, see Supplementary Figures 1A and 1B for a corresponding summaries of all SNPs.

Results were very similar for the analysis based on the T1DGC dataset adjusting for age, sex, region of birth and HLA risk category based on DQ genotypes (Supplementary Table 2, Figure 1).

ERBB3 gene and IA-2A

The association between the carrier status of the minor AA+AC genotype in the ERBB3 gene and T1D was modified by IA-2A (p for comparison=0.047), adjusting for age, sex, region of birth and HLA risk category (Table 4, Figure 1). The OR of association between ERBB3 (AA+AC) and IA-2A-positive T1D (OR=1.79 (1.02, 3.13); p=0.042) was higher compared to IA-2A-negative T1D (OR=1.28 (0.74, 2.21); p=0.378) (p for comparison=0.047) (Table 4, Figure 1).

Table 4

(Top) Estimated odds ratios (OR) and 95% confidence intervals (95% CI) of being diagnosed with T1D for subjects with the variant (minor) genotype compared to subjects with the major genotype and with low risk HLA haplotype (reference group), adjusting for age, sex and region of birth. The fold increase column summarizes the estimated OR (95% CI) of being diagnosed with T1D for carriers of the variant allele compared to non-carriers, for subjects in a given HLA risk category. The p-value = 0.047 for the model without interaction (top row), p-value = 0.056 for the model with interaction of ERBB3 and HLA risk category, and the p-value for trend = 0.176. (Bottom) The estimated OR and 95% CI of being diagnosed with autoantibody positive T1D (or autoantibody negative T1D) for subjects with the minor genotype compared to those with the major genotype, adjusted for age, sex, region of birth and HLA risk category.

HLA risk category	ERBB3 (CC)		ERBB3 (AA+AC)		Fold increase OR (95% CI)
	nP/nC	OR (95% CI)	nP/nC	OR (95% CI)
	101/81	(reference)	178/109	1.58 (1.01, 2.48)	1.58 (1.01, 2.48)
Low	11/35	(reference)	14/56	0.86 (0.34, 2.14)	0.86 (0.34, 2.14)
Neutral	22/26	3.19 (1.28, 7.96)	58/36	5.72 (2.54, 12.89)	1.79 (0.86, 3.72)
High	68/20	11.25 (4.77, 26.53)	106/17	22.65 (9.45, 54.30)	2.01 (0.98, 4.16)

Low HLA risk = {(DQB1*0602 or DQB1*0603), (DQB1*0301 but not DQB1*02, DQB1*0302, DQB1*0604)}

Neutral HLA risk = {(DQA1*0501-B1*0201 but not DQA1*0301-B1*0302), Other genotypes (not appearing in low or high categories)}

High HLA risk = {(DQA1*0501-B1*0201 and DQA1*0301-B1*0302), (DQA1*0301-B1*0302 but not DQA1*0501-B1*0201)}

Similar observation was made based on analysis using HLA risk categories based on DQ genotypes (Supplementary Table 3, Figure 1).

Gene 26471, IAA and IA-2A

The OR of association between the common GG genotype in the 26471 gene and T1D was OR=2.33 (1.44, 3.77); p=0.001 (Table 5, Figure 1).

Table 5

(Top) Estimated odds ratios (OR) and 95% confidence intervals (95% CI) of being diagnosed with T1D for subjects with the variant (major) genotype compared to subjects with the minor genotype and with low risk HLA haplotype (reference group), adjusting for age, sex and region of birth. The fold increase column summarizes the estimated OR (95% CI) of being diagnosed with T1D for carriers of the variant allele compared to non-carriers, for subjects in a given HLA risk category. The p-value = 0.001 for the model without interaction (top row), p-value = 0.002 for the model with interaction of 26471 and HLA risk category, and the p-value for trend = 0.636. (Bottom) Estimated OR and 95% CI of being diagnosed with autoantibody positive T1D (or autoantibody negative T1D) for subjects with variant (major) genotype compared to those with the minor genotype, adjusted for age, sex, region of birth and HLA risk category.

HLA risk category	26471 (AA+AG)		26471 (GG)		Fold increase OR (95% CI)
	nP/nC	OR (95% CI)	nP/nC	OR (95% CI)
	171/147	(reference)	108/42	2.33 (1.44, 3.77)	2.33 (1.44, 3.77)
Low	17/74	(reference)	9/17	2.29 (0.86, 6.11)	2.29 (0.86, 6.11)
Neutral	42/45	4.39 (2.2, 8.79)	41/16	12.55 (5.63, 27.99)	2.86 (1.38, 5.91)
High	112/28	18.22 (9.13, 36.38)	58/9	32.22 (13.07, 79.45)	1.77 (0.77, 4.04)

Low HLA risk = {(DQB1*0602 or DQB1*0603), (DQB1*0301 but not DQB1*02, DQB1*0302, DQB1*0604)}

Neutral HLA risk = {(DQA1*0501-B1*0201 but not DQA1*0301-B1*0302), Other genotypes (not appearing in low or high categories)}

High HLA risk = {(DQA1*0501-B1*0201 and DQA1*0301-B1*0302), (DQA1*0301-B1*0302 but not DQA1*0501-B1*0201)}

A strong relationship between gene 26471 and T1D was observed in the analysis using the HLA risk categories based on DQ genotypes. The overall OR of being diagnosed with T1D for those with the GG genotype compared to subjects with the AA+AG genotype was 2.79 (1.73, 4.50); p=2.46e-5) (Supplementary Table 4, Figure 1).

Furthermore, the results suggest that the association was stronger between IAA-negative T1D and 26471 gene (OR=2.74 (1.6, 4.66); p=2.17e-4) compared to IAA-positive T1D and the 26471 gene (OR=1.51 (0.78, 2.92); p=0.223) (p for comparison=0.049). Similarly for IA-2A, where the association appeared stronger between IA-2A-negative T1D and 26471 gene (OR=2.97 (1.66, 5.30), p=2.45e-4) compared to IA-2A-positive T1D and the 26471 gene (OR=1.47 (0.81, 2.69), p=0.205) (p for comparison=0.052) (Table 5, Figure 1). These results were supported by the analysis using the HLA risk categories based on DQ genotypes (Supplementary Table 4, Figure 1).

Major insulin gene AA genotype and IAA

The OR of association between the common AA genotype in the INS gene and T1D was significant (OR=2.00 (1.28, 3.11); p=0.002). In addition, the OR of association between the insulin gene and IAA-positive T1D patients (OR=4.16 (2.11, 8.21); p=4.03e-5) was higher compared to that for IAA-negative T1D patients (OR=1.6 (0.98, 2.63)) (p=0.060) and this difference in association was nominally statistically significant (p for comparison=0.008) (Table 6, Figure 1).

Table 6

(Top) Estimated odds ratios (OR) and 95% confidence intervals (95% CI) of being diagnosed with T1D for subjects with the (major) genotype compared to subjects with the minor genotype and with low risk HLA haplotype (reference group), adjusting for age, sex and region of birth. The fold increase column summarizes the estimated OR (95% CI) of being diagnosed with T1D for carriers of the variant allele compared to non-carriers, for subjects in a given HLA risk category. The p-value = 0.002 for the model without interaction (top row), p-value = 0.009 for the model with interaction of INS and HLA risk category, and the p-value for trend = 0.861. (Bottom) Estimated OR and 95% CI of being diagnosed with autoantibody positive T1D (or autoantibody negative T1D) for subjects with variant (major) genotype compared to those with the minor genotype, adjusted for age, sex, region of birth and HLA risk category.

HLA risk category	INS (AT+TT)		INS (AA)		Fold increase OR (95% CI)
	nP/nC	OR (95% CI)	nP/nC	OR (95% CI)
	82/93	(reference)	195/96	2.00 (1.28, 3.11)	2.00 (1.28, 3.11)
Low	8/49	(reference)	17/40	2.45 (0.94, 6.37)	2.45 (0.94, 6.37)
Neutral	27/25	6.76 (2.64, 17.31)	56/36	10.03 (4.18, 24.02)	1.48 (0.74, 2.97)
High	47/19	15.74 (6.13, 40.47)	122/20	38.36 (15.51, 94.85)	2.44 (1.18, 5.01)

Low HLA risk = {(DQB1*0602 or DQB1*0603), (DQB1*0301 but not DQB1*02, DQB1*0302, DQB1*0604)}

Neutral HLA risk = {(DQA1*0501-B1*0201 but not DQA1*0301-B1*0302), Other genotypes (not appearing in low or high categories)}

High HLA risk = {(DQA1*0501-B1*0201 and DQA1*0301-B1*0302), (DQA1*0301-B1*0302 but not DQA1*0501-B1*0201)}

The results were very similar in the analysis using the HLA risk categories based on DQ genotypes (Supplementary Table 5, Figure 1).

Major IL2RA AA genotype and IAA

The OR of association between the common AA genotype in the IL2RA gene and T1D was significant (OR=1.90 (1.09, 3.33); p=0.025). In addition, the OR of association between IL2RA and IAA-positive T1D patients (OR=3.91 (1.65, 9.28); p=0.002) was higher compared to that for IAA-negative T1D patients (OR=1.53 (0.83, 2.80); p=0.171) and this difference in association was nominally statistically significant (p for comparison=0.028) (Table 7, Figure 1).

Table 7

(Top) Estimated odds ratios (OR) and 95% confidence intervals (95% CI) of being diagnosed with T1D for subjects with the variant (major) genotype compared to subjects with the minor genotype and with low risk HLA haplotype (reference group), adjusting for age, sex and region of birth. The fold increase column summarizes the estimated OR (95% CI) of being diagnosed with T1D for carriers of the variant allele compared to non-carriers, for subjects in a given HLA risk category. The p-value = 0.025 for the model without interaction (top row), p-value = 0.042 for the model with interaction of IL2RA and HLA risk category, and the p-value for trend = 0.300. (Bottom) Estimated OR and 95% CI of being diagnosed with autoantibody positive T1D (or autoantibody negative T1D) for subjects with variant (major) genotype compared to those with the minor genotype, adjusted for age, sex, region of birth and HLA risk category.

HLA risk category	IL2RA (AG+GG)		IL2RA (AA)		Fold increase OR (95% CI)
	nP/nC	OR (95% CI)	nP/nC	OR (95% CI)
	41/39	(reference)	222/149	1.9 (1.09, 3.33)	1.9 (1.09, 3.33)
Low	3/13	(reference)	21/78	1.27 (0.32, 4.94)	1.27 (0.32, 4.94)
Neutral	13/14	4.63 (1.05, 20.46)	65/46	7.21 (1.91, 27.26)	1.56 (0.66, 3.68)
High	25/12	10.46 (2.43, 45.00)	136/25	27.48 (7.14, 105.81)	2.63 (1.15, 6.00)

Low HLA risk = {(DQB1*0602 or DQB1*0603), (DQB1*0301 but not DQB1*02, DQB1*0302, DQB1*0604)}

Neutral HLA risk = {(DQA1*0501-B1*0201 but not DQA1*0301-B1*0302), Other genotypes (not appearing in low or high categories)}

High HLA risk = {(DQA1*0501-B1*0201 and DQA1*0301-B1*0302), (DQA1*0301-B1*0302 but not DQA1*0501-B1*0201)}

The results were very similar in the analysis using the HLA risk categories based on DQ genotypes (Supplementary Table 6, Figure 1).

None of the other non-HLA genetic polymorphisms reported by T1DGC as nominally statistically significantly associated with T1D showed a relationship with T1D, T1D stratified by HLA-DQ or islet autoantibody-specific T1D in our dataset (Supplementary Tables 7A-13B).

Discussion

Supporting evidence for at least 24% of the non-HLA genes previously identified in the T1DGC [15, 16] to be associated with T1D was found in the present study of 305 T1D patients and 203 controls. The rank order of the genes associated with T1D was consistent with recent reports based on T1DGC data [16, 18]. More importantly, the second aim to test if the interrelation between low risk HLA-DQ and PTPN22, and that between PTPN22 and GADA-specific T1D [14] was successfully verified in this smaller dataset. Also, the observations that the major allele (A) of INS interacts with high risk HLA, and that it is associated with IAA-positive T1D, are consistent with previous findings [21, 33, 36, 37]. Additional major findings from our analyses were that: 1) ERBB3 increased the risk of IA-2A-positive T1D; 2) the risk carrying G allele of 26471 was associated with IAA-negative and IA-2A-negative T1D; 3) the risk carrying major allele A of IL2RA was associated with IAA-positive T1D and appears to double the risk of T1D among subjects in the high risk HLA group. These interactions at the time of the clinical diagnosis of T1D support the view that the T1D pathogenesis is highly heterogeneous. Thus, the statistical approach used in the previous [14] and in the present study should prove useful when larger numbers of T1D patients are studied.

In the analyses stratified by the HLA risk category (Table 2, top panels of Tables 3-7) we show not only the fold increase for each HLA risk category, but also the OR of the combined effect of HLA risk category and the SNP compared to the reference (low risk HLA and non-risk allele of a given SNP). Results presented that way clearly show the importance of HLA in the disease process, and that a fold increase of, for example, two attributed to a given SNP in a given HLA risk category has more severe consequences for patients in the high risk compared to low risk HLA. It will be equally important to apply the present approach to studying progression to T1D in children who have developed islet autoimmunity, defined as having developed persistent islet autoantibodies, which is the first sign of a beta-cell destructive process that after months to years may result in a clinical diagnosis of diabetes.

The PTPN22 T variant leads to a substitution of an arginine for a tryptophan residue at position 620 in the PTPN22 protein [38]. The mutation is thought to allow T cells to remain functionally active for a longer period of time. The OR of association between PTPN22T and GADA-positive T1D declined with increasing HLA risk while there was no such change in the GADA-negative T1D patients [14]. Our finding that the association between PTPN22T and T1D was modified by GADA and HLA, whether risk categories or actual genotypes were used to stratify subjects by risk of T1D, suggests that it will be important to explore antigen-presentation of GAD65 on HLA heterodimers in relation to T cell function. At present we do not have enough information to determine whether our results are due to a triggering mechanism of islet autoimmunity (step one), progression to clinical diagnosis of diabetes, or both. Taken together, our results suggest that the altered T-cell function conferred by the Arg620Trp substitution was exclusively associated with GADA-positive T1D. This observation should be taken into account both when children genetically predisposed to T1D are followed until the development of the first islet autoantibody be it IAA, GADA, IA-2A or, perhaps, ZnT8A [8, 39]. Hence, the PTPN22 polymorphism appears to be more important to the risk for seroconversion to GADA than IAA.

The ERBB3 gene is coding for a tyrosine kinase expressed as a cell surface growth factor receptor widely expressed in many developing mammalian tissues, including in the intestinal tract [40]. The risk of being diagnosed with IA-2A-positive T1D was increased in individuals with the minor (C) allele of ERBB3 gene. Longitudinal studies have suggested that IA-2A is rarely the first autoantibody to appear at the time of seroconversion to islet autoimmunity [8, 39, 4142, 43]. Screening for IA-2A and ZnT8A among parents and siblings was found to be a more cost-effective way to predict clinical onset than employing the entire battery of islet autoantibody tests [43]. We speculate that the minor ERBB3 C allele may be a useful marker for this type of IA-2A related progression to clinical onset.

A more complex relationship was observed between for the rare allele of the yet unknown gene 26471and IAA-negative T1D as well as possibly IA-2A-negative T1D. These interactions provide additional evidence suggesting that the T1D disease pathogenesis is highly heterogeneous.

The Swedish 1985-1989 dataset has been previously analyzed to evaluate the relationship between the HLA and islet autoantibodies [21], the GIMAP5 gene and IA-2 autoantibodies [44], the interaction between HLA-DQ, GAD65 autoantibodies and PTPN22 [14] and identification of ITPR3 gene as a possible non-HLA gene increasing the risk of diabetes [45]. Although the number of subjects was reduced due to prior extensive analyses, especially the use of larger DNA quantities in the pre-PCR era [46, 47], it was still possible to detect hitherto not previously recognized interactions: ERBB3 and IA-2A, 26471 and IAA, and possibly IA-2A, as well as IL2RA and IAA.

The samples which did not have an adequate amount of DNA required by the T1DGC seemed to represent a random sample compared to our most recent study of this 1985-1989 Swedish population-based case-control study [14]. Our sample size was large enough to support our previous findings stratifying on HLA risk group, as well as in relation to autoantibody-specific T1D. However, it was too small for us to try to consider the association between autoantibody-specific T1D and T1DGC identified SNPs stratifying by the HLA risk group, as we were able to do in our previous analysis [14]. In our analysis we were able to find supporting evidence for 24% of the SNPs identified by the T1DGC as nominally statistically significantly associated with T1D. Our inability to find evidence of an association between T1D and the remaining SNPs could be the result of a lack of power to detect smaller effect sizes detected in the T1DGC studies. Alternatively, the lack of replication for some SNPs could be due to the differences in the Swedish population studied by us and the multinational and more heterogeneous population studied by the T1DGC. It would be interesting to determine which T1DGC SNPs could be replicated in similar datasets from other countries. A larger Swedish dataset would allow us to test not only whether other SNPs could be replicated, but also whether the association of any of ERBB3, gene 26471, INS or IL2RA SNPs and islet autoantibodies are modified by the HLA risk group, as we previously reported for PTPN22, GAD65 autoantibody-positive T1D and low risk HLA [14].

In future studies with larger cohorts, such as the ongoing BDD study of consecutively diagnosed 1-18 year old patients now reaching more than 4000 patients [9, 28], should prove useful to further dissect the possible interaction between the T1DGC identified SNPs, HLA-DQ and islet autoantibodies. Also, further investigations are needed to determine if the T1DGC identified SNPs are primarily associated with the risk for islet autoantibodies in individuals with certain HLA-DQ rather than with the clinical onset of T1D. In the future, when selecting individuals for primary prevention trials, it will be important to understand whether the T1DGC identified SNPs and HLA-DQ need to be considered as risk factors of islet autoimmunity. Furthermore, could it be possible that for some non-HLA genes there is no association with T1D, and that such genes are risk factors for progression to T1D only in those subjects already positive for islet autoantibodies?

We conclude that in the present case-control study we found supporting evidence for the association between T1D and about 24% of the SNPs reported by the T1DGC [16]. We also found supporting evidence that the minor PTPN22 allele contributed to an increased T1D risk primarily in patients with neutral risk HLA-DQ. In addition, we report three novel findings. First, the ERBB3 gene and T1D association was modified by IA-2A. Second, the association between T1D and the 26471 gene was stronger both in IAA-negative and IA-2A-negative T1D patients. Third, the common INS genotype was associated with T1D patients positive for IAA but its association with T1D did not appear to be modified by HLA-DQ. Fourth, the association between IAA-positive T1D and IL2RA was nominally statistically significantly stronger than that between IAA-negative T1D and IL2RA, but we did not detect an interaction with HLA risk categories based on DQ genotypes. The SNPs evaluated in this study are likely to be of most value in studying the Swedish population, as the frequencies of genetic polymorphisms vary among populations. However, our analysis and results demonstrating that non-HLA genes affect T1D risk differently depending on the presence of islet autoantibodies and HLA risk categories are likely to be universal. Thus, these interactions should be studied in other populations and taken into account when subjects are randomized into both primary and secondary prevention trials, especially when autoantigens will be used to halt the progression to clinical onset of T1D.

Methods

Subjects

The present study was carried out in a subset of subjects from a matched case-control study of Swedish T1D patients (n=1006) and controls (n=702) matched based on age, gender and region of birth, as previously described in detail [14, 21, 47-49]. Due to prior extensive analyses of the original set of DNA samples [14, 21, 44, 45], the number of available DNA samples to fulfill the requirement of the T1DGC for SNP analysis on the Illumina typing system had diminished (Supplementary Table 14). Contingent on available DNA, in the present study it was possible to fully genotype 305 T1D patients and 203 controls. The National Ethics Board approved the study and all patients and controls gave informed consent.

In order to make maximum use of available data, all subjects with sufficient amount of DNA for genetic typing were used in an unmatched analysis (n=508), similarly to the previously reported dataset (n=1240) where PTPN22 was the SNP of interest [14] (Supplementary Table 14). To verify the appropriateness of this approach, we compared the results of conditional logistic regression analysis of the 148 subjects in matched sets with results of ordinary logistic regression analysis with adjustment for age, sex and region of birth. Since the regression coefficients were very similar (data not shown), we decided to account for the matching using adjustment rather than to perform a matched analysis.

Based on the HLA classification used in the DiPiS study[50], we defined three HLA risk categories (high, neutral and low) to account for the confounding effects of HLA [14]. The HLA risk categories (referred to as the original HLA risk classification in the text) were defined as follows: Low = {(DQB1*0602 or DQB1*0603), (DQB1*0301 but not DQB1*02, DQB1*0302, DQB1*0604)}, Neutral = {(DQA1*0501-B1*0201 but not DQA1*0301-B1*0302), Other genotypes (not appearing in low or high categories)}, High = {(DQA1*0501-B1*0201 and DQA1*0301-B1*0302), (DQA1*0301-B1*0302 but not DQA1*0501-B1*0201)} [14].

We also introduce a second HLA risk classification based on DQ genotypes [28, 51], where DQ2/8 and DQ8/8 are defined as high risk, DQ2/X as neutral, and all other DQ genotypes are defined as low risk.

We defined autoantibody positive and negative subjects based on GADA, IAA, IA-2A and ICA using previously defined cutoffs [21].

The main analysis and results described in this paper are based on 508 subjects, 305 T1D patients and 203 controls (Supplementary Table 14), referred to as the Swedish T1DGC (S-T1DGC) dataset.

HLA typing

HLA typing of DQA1 and DQB1 was performed by PCR amplification of the second exon of the genes followed by dot blot hybridizations of sequence-specific oligo probes and by restriction fragment-length polymorphism using DR- and DQ-based probes to establish haplotype [21, 46, 52]. HLA-DQ types were available on all 508 subjects in the S-T1DGC dataset.

Islet autoantibody analyses

Data was available on autoantibodies against insulin, GAD65 and IA-2. All serum and plasma samples have been expended in prior studies and it has therefore not been possible to determine the presence of the recently described ZnT8 transporter autoantibodies [53, 54].

Insulin autoantibodies (IAA)

Insulin autoantibodies (IAA) as well as antibodies to insulin were measured by radiobinding assay using acid–charcoal extraction and cold insulin displacement as described [55] IAA-Δ percent results were available on 196 of the 203 patients (96.6%) and 299 of the 305 control subjects (98.0%) in the S-T1DGC dataset.

GAD65 autoantibodies (GADA)

Autoantibodies to radiolabeled human Mr 65 000 GAD65 (GADA) were quantified by fluid-phase immunoprecipitation assay using fluorographic densitometry as described [56, 57] except that for 15- to 35-year-olds, GAD65 were radiolabeled by coupled in vitro transcription and translation, as described [58, 59]. The two assay formats correlate well, and autoantibody levels are expressed as GADA index [59], using the WHO–Juvenile Diabetes Foundation (JDF) standard [60] as the positive control. GADA measurements were available on 200 of the 203 patients (98.5%) and 297 of the 305 control subjects (97.4%) in the ST1DGC dataset.

Islet antigen-2 (IA-2) autoantibodies (IA-2A)

Antibodies to IA-2 [61, 62] were measured by radiobinding immunoassay [63]. In vitro translation with 35S-methionine yielded a 46-kDa IA-2 polypeptide highly precipitable by diabetes sera. Radiobinding assays used protein A-Sepharose to separate antibody bound from free IA-2, and autoantibody levels are expressed as an IA-2A index [59], using the WHO–JDF standard as a positive control [60]. IA-2A measurements were available on 198 of the 203 patients (97.5%) and 277 of the 305 control subjects (90.8%) in the ST1DGC dataset.

Islet cell (cytoplasmic) antibodies (ICA)

ICA were determined by indirect two-color immunofluorescence using blood group O frozen human pancreas [64]. The same pancreas was used throughout the study. Two independent observers evaluated coded slides. Samples were titered in doubling dilutions to determine end points for conversion into JDF units [65] as described [64]. ICA of 0 JDF units corresponded to no detectable antibody. ICA measurements were available on 203 of the 203 patients (100%) and 304 of the 305 control subjects (99.7%) in the S-T1DGC dataset.

Statistical analyses

The analyses performed for the S-T1DGC SNPs were based on those performed in [14] for the rs2476601 SNP in PTPN22. Using the dbSNP database (http://www.ncbi.nlm.nih.gov/SNP/) we identified the common allele ‘M’ (referred to as ancestral or major) and the rare (minor) allele ‘m’ for each SNP. Regardless of whether the risk carrying allele was common or rare, we combined the ‘mm’ genotype with the ‘Mm’ genotype due to low numbers of subjects in the ‘mm’ group for most SNPs. In future studies it will be important to extend the SNP analysis to the three SNP genotypes. It cannot be excluded that mode of inheritance of a heterozygous or homozygous predisposing allele may affect the combined interaction between a SNP, HLA-DQ and islet autoantibodies.

The variables in our dataset were coded as follows. For each SNP, if ‘M’ was the risk carrying allele, ‘MM’ was coded as 1 and ‘Mm+mm’ as 0; if ‘m’ was the risk carrying allele, ‘mm+Mm’ was coded as 1 and ‘MM’ as 0; age as a continuous variable, region as unordered categorical, and the autoantibody status for GADA, IAA, IA-2A, and ICA as binary. HLA haplotypes were grouped into low-, neutral-, and high-risk HLA, and these HLA-risk categories were encoded in two ways, as an ordered categorical variable (low, neutral, and high), and as a linear effect (−1, 0, 1). The outcome, or the case–control (T1Dcc) status, was coded as a binary variable (case = 1, control = 0).

In the analysis, we considered all T1D cases, as well as what we refer to as “an autoantibody-positive (or autoantibody-negative) T1D”, which were subgroups of T1D patients who were positive (or negative) for a given autoantibody, regardless of their status with respect to other autoantibodies. In this way, we created eight additional data subsets. In each of these, we included all of the controls (regardless of their autoantibody status), and one of the following: autoantibody-positive T1D cases, autoantibody-negative T1D cases, in which autoantibody was one of GADA, IAA, IA-2A, or ICA. Hence, autoantibody status is used to divide T1D cases into subsets, as is commonly performed in RA and other autoimmune diseases [14].

For SNPs with ‘M’ being the risk carrying allele, all measures of association refer to the odds ratios (OR) of association for subjects with the ‘MM’ genotype and T1D compared to subjects with the ‘Mm+mm’ genotype. For SNPs with ‘m’ being the risk carrying allele, all measures of association refer to the ORs of association in subjects with the ‘mm+Mm’ genotype and T1D compared to subjects with the ‘MM’ genotype.

For each SNP in the T1DGC dataset we performed the following analysis: Step 1 - We estimated the OR of association between the risk carrying allele of a given SNP (rSNP) and T1D by fitting a logistic regression model with T1D as the outcome, with the SNP as the main effect, and adjusting for age, sex, region (design variables), and the HLA risk category. Step 2 — To estimate the OR in each of the HLA risk categories, we fit a model similar to that in Step 1, but with an added interaction term between the HLA and the SNP. Step 3 — To determine whether the association between the rSNP and T1D was associated with the appearance of autoantibodies we fit the models described in Step 1 to the autoantibody-positive T1D dataset (note: all controls were included in all subsets of the data), similarly for autoantibody-negative T1D dataset, in which autoantibody was one of GADA, IAA, IA-2A, or ICA. Step 4 — To test whether the association between the rSNP and T1D was nominally statistically significant, we used the Wald’s test with one degree of freedom (d.f.); similarly for autoantibody-positive T1D and autoantibody-negative T1D, in which autoantibody was one of GADA, IAA, IA-2A, or ICA. Step 5 — To compare the estimated ORs of association between the rSNP and autoantibody-positive T1D and that between the rSNP and autoantibody-negative T1D, we used a Wald's test with d.f.=1 based on estimates of mean and variance from polytomous logistic regression[66] . It allowed us to model an outcome with three levels: controls, autoantibody-positive T1D cases, and autoantibody-negative T1D cases, thus accounting for the correlation between the analyses involving autoantibody-positive T1D and autoantibody-negative T1D, in which the correlation was due to the fact that the same set of controls was used in the pair of analyses for a given autoantibody. Step 6 — We considered whether there was a linear trend in the interaction effect of HLA-risk category on the association between the rSNP and T1D. To accomplish that, we fit the models described in Steps 2 and 3, adjusting for HLA risk (coded as a categorical variable as before), but now using the HLA coded as linear effect in the interaction term. We used a likelihood ratio test with d.f. = 1 to test the appropriateness of the linear HLA effect assumption. We evaluated the appropriateness of our assumptions using residual plots for the models we fitted.

All of these analyses were performed for all SNPs for both HLA risk classifications. Thus, there were two tests performed for each type of analysis for each SNP, due to adjustment for each of the two HLA risk classifications. We report the results of all analyses (in the manuscript and Supplementary information), with the p-values not adjusted for multiple comparisons, for that reason we refer to p-values as “nominally statistically significant”.

All statistical analyses were performed in R v.2.12.1 (http://www.r-project.org). Data and the R code used in the analyses are available upon request.