Seven newly identified loci for autoimmune thyroid disease.

Autoimmune thyroid disease (AITD), including Graves' disease (GD) and Hashimoto's thyroiditis (HT), is one of the most common of the immune-mediated diseases. To further investigate the genetic determinants of AITD, we conducted an association study using a custom-made single-nucleotide polymorphism (SNP) array, the ImmunoChip. The SNP array contains all known and genotype-able SNPs across 186 distinct susceptibility loci associated with one or more immune-mediated diseases. After stringent quality control, we analysed 103 875 common SNPs (minor allele frequency >0.05) in 2285 GD and 462 HT patients and 9364 controls. We found evidence for seven new AITD risk loci (P < 1.12 × 10(-6); a permutation test derived significance threshold), five at locations previously associated and two at locations awaiting confirmation, with other immune-mediated diseases.

INTRODUCTION

Genome-wide association (GWA) studies have had a dramatic impact on susceptibility locus discovery and in addition, highlighted and extended the previously observed (1,2) commonality of loci between immune-mediated diseases, which has ultimately resulted in the ImmunoChip project. As described previously (3), the ImmunoChip project is a collaboration between 12 immune-mediated disease groups [autoimmune thyroid disease (AITD), ankylosing spondylitis (AS), celiac disease (CeD), Crohn's disease (CD), IgA deficiency, multiple sclerosis (MS), primary biliary cirrhosis (PBC), psoriasis (PS), rheumatoid arthritis (RA), systemic lupus erythematosus (SLE), type 1 diabetes (T1D) and ulcerative colitis (UC)] and the Wellcome Trust Case Control Consortium (WTCCC). The main objective of the project was to provide cost-effective genotyping of all known and genotype-able single-nucleotide polymorphisms (SNPs) across confirmed and distinct susceptibility loci from the 12 disease groups (www.ImmunoBase.org). The result was the ImmunoChip, a custom Illumina 200 K Infinium high-density array containing 186 distinct susceptibility loci (3), mainly derived from GWA studies, which have been convincingly associated (required P < 5 × 10−8) with one or more immune-mediated diseases. In addition, the ImmunoChip contained wildcard SNPs provided by the 12 disease groups and SNPs for replication found in non-immunological disease GWA studies conducted by the WTCCC (www.wtccc.org.uk/ccc2/wtccc2_studies.shtml).

The AITDs are caused by an autoimmune attack in which antibodies directed against components of the thyroid gland lead to either thyroid hormone excess and clinical features of hyperthyroidism as seen in people with Graves' disease (GD) or thyroid gland damage, reduced thyroid hormone production and clinical features of hypothyroidism as seen in people with Hashimoto's thyroiditis (HT). AITD shows a strong female preponderance with estimated incidence rates of GD in females as 80 per 100 000 person-years and in males as 8 per 100 000 person-years, and of HT in females as 350 per 100 000 person-years and in males as 8 per 100 000 person-years (4).

Although many regions of the genome have been reported as AITD susceptibility loci, convincing evidence for association has been limited to four regions: PTPN22/1q13.2, CTLA4/2q33.2, HLA/6p21, and TSHR/14q31.1 (5). These susceptibility loci all fall within the antigen presentation/T-cell receptor signalling pathway. In addition, there is mounting evidence for a susceptibility locus in FCRL3/1q23.1 (5–7), SCGB3A2/5q32 (8,9) and IL2RA/10p15.1 (10,11), and a recent GWA and replication study in subjects recruited from the Chinese Han population identified two novel loci at 4p14 and RNASET2-FGFR1OP-CCR6/6q27 (5).

Here, we used the ImmunoChip to further investigate the genetic architecture of AITD in 2374 GD and 474 HT patients, and 9953 controls.

RESULTS

After consideration of the number of AITD patients available, we focused on susceptibility locus discovery using common SNPs [minor allele frequency (MAF) ≥0.05] outside of the HLA, which reduced the number of tests performed and consequently, allowed us to use a less stringent significance threshold. We adopted a permutation approach to correct the 0.05 level of significance for multiple testing and for correlation between SNP genotypes, and obtained a threshold P = 1.12 × 10−6 (see Materials and Methods). After stringent quality control (see Materials and Methods), we analysed 103 875 common SNP genotypes in 2285 GD and 462 HT patients, and 9364 controls. As the vast majority of our AITD patients have been previously analysed in either discovery or confirmation studies of the reported AITD susceptibility loci, we initially used the AITD susceptibility loci on the ImmunoChip as positive controls and as expected, found evidence of association (P < 1.12 × 10−6 and good genotype signal intensity clusters; Supplementary Material, Table S1) for PTPN22 (P = 9.7 × 10−23), CTLA4 (P = 2.1 × 10−23), IL2RA (P = 2.7 × 10−7) and TSHR (P = 1.3 × 10−38 in GD patients/controls). We note that the AITD risk locus SCGB3A2 was not on the ImmunoChip, but had previously been associated in our AITD patients (P = 0.007; (8)). In addition, we found some evidence of association for FCRL3 (P = 1.1 × 10−5 in GD patients/controls; Supplementary Material, Table S1) and found further evidence of association for the recently reported RNASET2–FGFR1OP–CCR6/6q27 locus (5) (imm_6_167 338101 P = 1.6 × 10−7; r2 between imm_6_167338101 and rs9355610 (5) in controls = 0.72; Table 1; Supplementary Material, Fig. S1) in subjects with European ancestry. The recently reported 4p14 locus (5) was not on the ImmunoChip.

Table 1.

Single locus test evidence for the seven newly identified AITD non-HLA susceptibility loci and 6q27 on the ImmunoChip.

Karyotype band	Most disease-associated SNP	Located within the gene	MAF in controls	Immune-mediated diseases previously associated with regions (P < 5 × 10−8) (reported SNP; r2 with reported SNP in controls)	2285 GD patients and 9364 controls		462 HT patients and 9364 controls		Heterogeneity in disease association	2282 GD, 451 HT patients and 9364 controlsa
					OR (95% CI)	P-value	OR (95% CI)	P-value	P-value	OR (95% CI)	P-value
1p36.32	rs2843403 C>T	MMEL1 (32 kb upstream of TNFRSF1)	0.362	CeD (26) (rs3748816/TNFRSF14-MMEL1; r2 = 0.99)PBC (27) (rs3748816; r2 = 0.99)RA (28) (rs3890745/TNFRSF14-MMEL1; r2 = 0.90)	0.84 (0.79–0.90)	7.94 × 10−7	0.97 (0.85–1.12)	0.696	0.0593	0.86 (0.81–0.92)	5.63 × 10−6
2p25.1	rs1534422 A>G	216 kb upstream of TRIB2	0.455	Unconfirmed T1D locusb	1.16 (1.09–1.24)	4.69 × 10−6	1.24 (1.08–1.41)	1.64 × 10−3	0.394	1.17 (1.11–1.25)	1.76 × 10−7
3q27.3/3q28	rs13093110 C>T	LPP	0.452	CeD (26) (rs1464510; r2 = 0.98)	1.18 (1.10–1.26)	8.17 × 10−7	1.20 (1.05–1.37)	7.09 × 10−3	0.797	1.19 (1.12–1.26)	3.69 × 10−8
6q15	rs72928038 G>A	BACH2	0.177	CD (29) (rs1847472; r2 = 0.39)CeD (26) (rs10806425/BACH2-MAP3K7; r2 = 0.30)T1D (21) (rs11755527; r2 = 0.19)c	1.21 (1.12–1.32)	3.63 × 10−6	1.30 (1.11–1.53)	1.36 × 10−3	0.417	1.23 (1.14–1.32)	1.23 × 10−7
6q27	imm_6_167338101 A>C	FGFR1OP (118 kb upstream of CCR6)	0.408	AITD (5) (rs9355610; r2 = 0.72)CD (30) (rs2301436/CCR6; r2 = 0.61)RA (31) (rs3093023/CCR6; r2 = 0.25Vitiligo (32) (rs2236313; r2 = 0.22)	0.84 (0.79–0.90)	3.30 × 10−7	0.88 (0.76–1.00)	0.0564	0.575	0.85 (0.80–0.90)	1.64 × 10−7
11q21	rs4409785 T>C		0.173	Unconfirmed MS locusd	1.21 (1.11–1.31)	5.37 × 10−6	1.34 (1.14–1.57)	3.54 × 10−4	0.234	1.23 (1.14–1.33)	7.69 × 10−8
12q12	rs4768412 C>T	PRICKLE1	0.363	CD (30) (rs11175593/LRRK2-MUC19)d	1.19 (1.11–1.27)	3.30 × 10−7	1.00 (0.88–1.15)	0.949	0.0234e	1.16 (1.09–1.23)	3.56 × 10−6
16p11.2	rs57348955 G>A	83 kb upstream of ITGAM	0.396	SLE (33) (rs1143679/ITGAM; r2 = 0.056)T1D (13) (rs4988084/IL27)f	0.83 (0.77–0.89)	3.76 × 10−8	0.91 (0.80–1.05)	0.188	0.190	0.84 (0.79–0.89)	5.13 × 10−8

The reported SNPs have a P-value of <1.12 × 10−6 (see Materials and Methods) and were the most disease-associated SNP in the region with good genotype signal intensity plots (www.ImmunoBase.org). For each SNP, the most significant P-value is shown in bold. Note that the same controls were analysed in each analysis. AITD, autoimmune thyroid disease; MAF, minor allele frequency; HT, Hashimoto's thyroiditis; GD, Graves' disease; OR, Odds ratio for the minor allele; CI, confidence interval; CD, Crohn's disease; CeD, celiac disease; PBC, primary biliary cirrhosis; and SLE, systematic lupus erythematosus.

aWe excluded 14 second degree relatives or closer between disease groups (see Materials and Methods).

brs1534422 at 2p25.1 is an unconfirmed T1D risk locus (T1D P = 2.1 × 10−6; (13)).

cWe used CEU data from the 1000 Genomes Project to estimate r2 as the T1D SNP, rs11755527, was not available in the ImmunoChip data.

drs4409785 at 11q21 is an unconfirmed MS risk locus [MS P = 6.3 × 10−7; (14)].

eAssuming 13 [8 (Table 1) and 5 (Supplementary Material, Table S1)] independent tests, the adjusted P-value was 3.85 × 10−3 for the 0.05 level of significance based on the Bonferroni correction for multiple testing.

fNo available genotype data which included both SNPs.

We found evidence of association between AITD and seven genomic locations (P < 1.12 × 10−6; Table 1). Five of these new AITD risk loci were at locations previously associated with at least one other immune-mediated disease (Table 1; Supplementary Material, Figs S2–S6); three of which included the most AITD-associated SNPs being located within MMEL1/1p36.32, LPP/3q28 and BACH2/6q15 genes. However, only BACH2 has a known autoimmune function as a regulator of nucleic acid-triggered antiviral responses in human cells (12). The remaining two new AITD risk loci were unconfirmed risk loci: rs1534422 at 2p25.1 is an unconfirmed T1D risk locus (T1D P = 2.1 × 10−6; (13); Supplementary Material, Fig. S7); and rs4409785 at 11q21 is an unconfirmed MS risk locus (MS P = 6.3 × 10−7; (14); Supplementary Material, Fig. S8).

We tested for heterogeneity in disease association before combining the GD and HT patients (Table 1 and Supplementary Material, Table S1). After allowing for multiple testing, only the disease association of the TSHR/14q31.1 locus was significantly different between the two patient groups. This difference between GD and HT patients was expected as the TSHR represents the primary autoantigen in hyperthyroidism of GD (15). As we focused on susceptibility loci discovery using common SNPs (MAF > 0.05) outside of the HLA to reduce the number of tests performed, which was necessary given the number of GD and HT patients available in this study, we did not proceed to test whether the most disease-associated SNP alone was sufficient to model the association signal in AITD risk loci with dense SNP coverage.

We examined the overlap between AITD and the other 11 immune-mediated diseases that provided susceptibility loci to the ImmunoChip in turn, by testing whether the pattern of AITD association across immune-mediated disease regions differed according to whether those regions had previously been associated with each specific immune-mediated disease or not (Table 2; (16)). We restricted the analysis to non-HLA SNPs from densely mapped regions associated with at least one immune-mediated disease. We found that SNPs within RA-associated regions were significantly (P = 6.56 × 10−5) more likely to be associated with AITD than SNPs within regions not associated with RA. In addition, there was marginal evidence of an overlap between AITD and T1D-associated regions (P = 0.0327). In contrast, we found that SNPs within regions associated with AS, UC and PS showed significantly less evidence of association with AITD than SNPs within regions known to be associated with each disease (P = 3.27 × 10−4, 1.69 × 10−3 and 4.08 × 10−3, respectively). There was also marginal evidence of an overlap between AITD and regions not associated with CD (P = 0.0250). We note that an important caveat about interpreting these test results is that the susceptibility loci for each disease are incomplete and the extent of the incompleteness varies between diseases. This prevents us, for example, from drawing any conclusions that ‘AITD is less like AS than UC’. However, for the majority of the immune-mediated diseases listed in Table 2, we can assume that the common loci with the largest effects are known and that significant results here are true representations of AITD-disease overlap.

Table 2.

The overlap between AITD and the 11 other immune-mediated diseases that provided susceptibility loci to the ImmunoChip. A negative Z-value indicates that SNPs within disease-associated regions were more likely to be associated with AITD than SNPs within regions not known to be associated with the disease. A positive Z-value indicates that SNPs within regions not known to be associated with the disease were more likely to be associated with AITD than SNPs within associated regions.

Immune-mediated disease	Number of non-HLA loci	Number of SNPs within loci	Wilcoxon rank test
			Z	P
AS	8	2079	3.593	3.27 × 10−4
CeD	30	9866	0.147	0.883
CD	66	22 585	2.242	0.0250
IgA deficiency	1	122	−1.208	0.227
MS	26	7147	−0.277	0.782
PBC	1	209	−0.840	0.401
PS	25	6767	2.872	4.08 × 10−3
RA	25	7742	−3.992	6.56 × 10−5
SLE	35	9099	−0.726	0.468
T1D	41	12 658	−2.146	0.0327
UC	28	8901	3.140	1.69 × 10−3

DISCUSSION

We have found evidence for seven risk loci not previously associated with AITD and expect that even more AITD risk loci will be discovered as more AITD patients are genotyped on the ImmunoChip. In this study, as the ImmunoChip contains densely mapped immune-mediated disease loci representing only a small proportion, ∼1%, of the genome, we used a permutation approach to derive a significance level rather than adopting a genome-wide significance level. Importantly, regardless of which discovery significance level was adopted, confirmation of the new AITD risk loci will require independent replication evidence. We found some statistical evidence that AITD, outside of the HLA, has a shared genetic architecture with both RA and T1D. AITD often co-exists with other immune-mediated diseases and a recent study estimating the prevalence of these co-existing disorders found that RA was the most common, occurring in 3.15% of GD patients and 4.24% of HT patients studied; T1D occurred in 1.11% of GD patients and 1.01% of HT patients (17). We note that the prevalence study (17) included vast majority of our AITD patients, and that AITD patients were recruited from endocrine clinics and not rheumatology clinics.

The implicated candidate genes for the established and seven new AITD risk loci indicate the importance of T lymphocyte signalling (PTPN22), T-regulatory cell function (CTLA4 and IL2RA), lymphocyte trafficking (CCR6) and a newly identified possible link with the anti-viral immune response (BACH2). Nevertheless, other AITD-associated regions suggest candidate genes of uncertain or unknown functions or roles in autoimmunity of AITD, such as LPP or TRIB2. Two regions, 11q21 and 12q12, have no annotated protein-coding candidate genes, indicating how much of the aetiology of this and other immune-mediated diseases remain to be explained.

MATERIALS AND METHODS

Subjects

The British AITD cases consisted of 2374 (390 males and 1984 females) GD and 474 (69 males and 405 females) HT patients (18). The British control population consisted of 6894 subjects drawn from the British 1958 Birth Cohort (1958BC; http://wwww.cls.ioe.ac.uk/studies.asp?section=000100020003) and 3059 subjects drawn from the UK Blood Services Common Control Collection (UKBS-CC) (6,19). All subjects were of white European ancestry with written informed consent and Ethics Committee/Institutional Review Board approval.

SNP selection

The dense SNP map for 186 confirmed and distinct loci, mainly derived from GWA studies, which have been associated (required P < 5 × 10−8) with one or more autoimmune diseases, consisted of all known SNPs, and small insertion deletions, in the dbSNP database, in the 1000 Genomes project (February 2010 release) and in additional sequencing data provided by collaborators. The ImmunoChip contains 196 524 polymorphisms (718 insertion/deletions and 195 806 SNPs) (3,20).

Genotyping

Samples were genotyped using a custom Illumina 200K Infinium high-density array according to the manufacturer's protocol at the Wellcome Trust Sanger Institute (WTSI) in Hinxton, UK, and University of Virginia (UVA) in Charlottesville, USA. In WTSI, 2374 GD, 474 HT, 1478 1958BC and 3059 UKBS-CC samples were processed, and in UVA, 5416 1958BC samples were processed.

Sample quality control

Samples were excluded based on per-sample call rate (Supplementary Material, Figs S9–S11), outlying autosomal heterozygosity (Supplementary Material, Figs S9–S11), inconsistent recorded and genotype-inferred sex, non-European ancestry (21) (Supplementary Material, Figs S12–S14), duplication and being closely related to another sample in the study. Relatedness between the study participants was estimated by a identity by state (IBS) statistic, we defined duplicates, or monozygotic twins, as having an IBS >0.98 and second degree relatives or closer, as having an IBS >0.1875 (22); the sample with the lower call rate was excluded from further analyses. Finally, we checked control sample's identity using the 169 intentional duplicates between the WTSI and the UVA, and using in-house genotype data. We identified 256 samples, which had incorrect identifiers as a result of manifest misreading, plate rotation and chip-related issues, and excluded 95 samples from unknown subjects (i.e. could not be assigned to identifiers). In addition, we found ten additional duplicates between the genotyping centers. All potential sample identity errors were discussed with the relevant genotyping centers by Neil M. Walker. After sample quality control, 2285 GD, 462 HT, 6541 1958BC and 2823 UKBS-CC samples were available for analysis. When we analysed GD and HT patients together, we excluded 14 second degree relatives or closer between the disease groups, which resulted in 2282 GD and 451 HT samples.

SNP quality control

SNP genotypes were called using Illumina GenomeStudio GenTrain 2.0 algorithm and genotypes with a gencall score (a quality metric for the distance a sample is from the center of the nearest genotype cluster) <0.15 were not assigned genotypes. SNPs were excluded if the MAF fell <5% in controls, if they deviated from Hardy–Weinberg equilibrium (HWE; z-score >5) in controls or if the per-SNP call rate fell <95% in cases or controls. We note that SNP quality control was performed in WTSI and UVA-processed controls separately. In addition, for SNPs from the pseudoautosomal regions on chromosome X, we tested for deviation from HWE in male and female controls separately after discovering SNP genotype frequency differences between sexes in patients and controls (e.g. rs306891). We excluded SNPs from the HLA (NCBI build 36 coordinates chr6:29690000-33498585).

Previously, despite differences in subject's age and in DNA processing between 1958BC and UKBS-CC control subjects, few differences in SNP allele frequencies between 1958BC and UKBS-CC controls were observed (19), which justified the combination of these two control groups into a single control group. To assess the possible differential genotyping errors between the WTSI and UVA genotyping centers, we tested for differences in SNP allele frequencies in controls between genotyping centers; the test was stratified by geographical region of Great Britain (23). The quantile–quantile plot for the WTSI control versus UVA control 1 degree-of-freedom (1-df) association tests (Supplementary Material, Fig. S15) showed a slight inflation of the test statistic [inflation factor (λ) was 1.017; as the inflation factor scales with sample size, a more informative and comparable measure of inflation is the inflation factor for an equivalent study of 1000 cases and 1000 controls, λ1000 (24), which was 1.004 for the WTSI control/UVA control analysis]. As the ImmunoChip was designed to include all known and genotype-able SNPs across 186 distinct immune-mediated disease susceptibility loci, there was inflation in the test statistic for the AITD patient versus control association tests: HT patient/control λ1000 = 1.06 (λ = 1.05); GD patient/control λ1000 = 1.08 (λ = 1.28); and HT and GD patient/control λ1000 = 1.08 (λ = 1.33). We note that λ was based on non-HLA common SNPs from regions not previously associated with AITD (excluded regions reported in Supplementary Material, Table S1) that passed quality control. We also estimated the inflation factor based on the WTCCC replication SNPs for reading and maths ability, and psychosis endophenotypes: HT patient/control λ1000 = 1.02 (λ = 1.02); GD patient/control λ1000 = 1.03 (λ = 1.11); and HT and GD patient/control λ1000 = 1.04 (λ = 1.17). We note that 2605 WTCCC replication SNPs passed quality control and had a MAF >0.05, and that the observed inflation of the test statistic was a result of the replication SNPs overlapping some of the immune-mediated disease loci, for example, the T1D loci COBL/7p12.1 and GLIS3/9p24.2. As the inflation factor λ1000 of the test statistics from the WTSI control/UVA control analysis was just over 1 and we found a similar inflation of the test statistics for GD patient/controls and for HT patient/controls at 6–8% in the main analysis and at 2–3% in the subset of WTCCC replication SNPs, population structure is unlikely to a major source of inflation. After SNP quality control, 103 875 common SNPs were available for analysis.

Statistical analysis

All statistical analyses were performed in R (http://www.r-project.org) using the snpStats package available from the Bioconductor project (http://www.bioconductor.org). After applying sample and SNP quality control exclusions, we analysed the case/control data using the 1-df Cochran-Armitage trend test. SNPs were modelled using a numerical indicator variable coded 0, 1 and 2, representing the occurrences of the minor allele and assuming multiplicative allelic effects. We also tested the disease-associated SNPs reported in Table 1 for differences between the multiplicative allelic effects model and the full genotype model, and no significant differences were found. We analysed SNPs on the X-chromosome using the method proposed by Clayton (25).

We initially analysed the GD and HT patients separately, but against the same controls. Before combining the GD and HT patients into a single patient group, as we had analysed GD and HT patients against the same controls, we tested for disease association heterogeneity between the two AITDs by testing for differences in SNP allele frequencies between GD and HT patient groups. In the absence of disease association heterogeneity, we proceeded to combine the GD and HT patients and re-analysed.

As GD and HT patients, unlike the controls, were not recruited from all 12 geographical regions of Great Britain (seven and six regions, respectively), we did not stratify the analysis by geographical region to allow for population structure. We note that previously, little difference was reported between analyses with or without regional stratification (19). However, to check that we were gaining power with the additional controls rather than false-positive findings through ignoring the population structure, we tested for population heterogeneity in SNP genotype frequencies across the regions to identify SNPs that differentiate between the regions (Supplementary Material, Table S2). None of the seven new or previously reported AITD risk loci showed large allele frequency differences between controls from the 12 regions.

To adjust P-values for the multiple tests and to allow for the correlation between SNP genotypes, we adopted a permutation test approach. We performed 100 000 permutation tests of the common SNPs on the ImmunoChip, permuting the disease status label for each test, to provide the distribution of the test statistic under the null hypothesis of no association between SNP genotypes and disease status. The adjusted P-value for 0.05 was 1.12 × 10−6.

To examine the overlap between AITD and the other 11 immune-mediated diseases that provided susceptibility loci to the ImmunoChip, we tested whether the distribution of AITD test statistics varied between SNPs within regions associated with another immune-mediated disease and within regions not associated with that disease. We adopted a gene set enrichment approach based on the Wilcoxon rank test as implemented in the R package wgsea (Chris Wallace, www-gene.cimr.cam.ac.uk/staff/wallace) (16). SNPs within the non-HLA immune-mediated disease loci were divided into two groups for each immune-mediated disease: (i) SNPs within regions associated with disease at the time the ImmunoChip was designed; and (ii) SNPs within regions not associated with that disease. We then tested the null hypothesis that the distribution of P-values describing the evidence for association of each SNP with AITD was the same in the two groups using a Wilcoxon rank test with a permutation-based variance estimate to allow for the considerable between SNP correlation (16). We used 10 000 permutations to estimate the variance of the Wilcoxon rank test.

SUPPLEMENTARY MATERIAL

Supplementary Material is available at .

FUNDING

The research leading to these results has received funding from the European Union's 7th Framework Programme (FP7/2007–2013) under grant agreement No. 241447 [Natural immunomodulators as novel immunotherapies for type 1 diabetes (NAIMIT)]. Funding to pay the Open Access publication charges for this article was provided by the Wellcome Trust.