A Genome-Wide Meta-Analysis of Six Type 1 Diabetes Cohorts Identifies Multiple Associated Loci

1Diabetes impacts approximately 200 million people worldwide, with microvascular and cardiovascular disease being the primary complications. Approximately 10% of cases are type 1 diabetes (T1D) sufferers, with ∼3% increase in the incidence of T1D globally per year. It is expected that the incidence is 40% higher in 2010 than in 1998.

2T1D is a clear example of a complex trait that results from the interplay between environmental and genetic factors. There are many lines of evidence that there is a strong genetic component to T1D, primarily due to the fact that T1D has high concordance among monozygotic twins and runs strongly in families, together with a high sibling risk.

3Prior to the era of GWAS, only five loci had been fully established to be associated with T1D. However, the majority of the other reported associations in the pre-GWAS era – remain highly doubtful, where an initial report of association does not hold up in subsequent replication attempts by other investigative groups. This previous hazy picture of the genetics of T1D can be put down to the use of the only methodologies that were available at the time and which were much more limited than GWAS i.e. the candidate gene approach (where genomic regions were studied based on biological reasoning) and family-based linkage methodologies. Inconsistent findings can also be attributed to small sample sizes i.e. when power is low the false discovery rate tends to be high; GWAS per se has not improved consistency, rather it has leveraged large, well powered sample sizes combined with sound statistical analyses.

4It has been long established that approximately half of the genetic risk for T1D is conferred by the genomic region harboring the HLA class II genes (primarily HLA-DRB1, -DQA1 and -DQB1 genes), which encode the highly polymorphic antigen-presenting proteins. Other established loci prior to the application of GWAS are the genes encoding insulin (INS) -, cytotoxic T-lymphocyte-associated protein 4 (CTLA4) –, protein tyrosine phosphatase, non-receptor type 22 (PTPN22) gene,, interleukin 2 receptor alpha (IL2RA) – and ubiquitin-associated and SH3 domain-containing protein A (UBASH3A).

5The application of genome wide association studies (GWAS) has robustly revealed dozens of genetic contributors to T1D –, the results of which have largely been independently replicated –. The most recently reported meta-analysis of this trait identified in excess of forty loci, including 18 novel regions plus confirmation of a number of loci uncovered through cross-disease comparisons –. As such, the risks conferred by these additional loci are relatively modest compared to the ‘low-hanging fruit’ described in the first studies and could only be ultimately uncovered when larger sample sizes were utilized.

6We sought to expand further on this mode of analysis by combining our cohort with all publically released genome wide SNP datasets to identify additional loci contributing to the etiology of T1D. Unfortunately, there is a relative paucity of control genotype data in these publically available sources. To circumvent this problem, we combined individual level data from each available cohort and we then compared the cases with controls from two sources. We next separated all the individual level data into two groups, characterized by the type of genotyping platform that was used to genotype the samples, which would later be recombined using inverse-variance meta-analysis. The 6,523 cases genotyped on an Illumina BeadChip included subjects from McGill University, The Children's Hospital of Philadelphia (CHOP), The Diabetes Control and Complications Trial – Epidemiology of Diabetes Interventions and Complications (DCCT-EDIC) cohort, and the Type 1 Diabetes Genetics Consortium (T1DGC), which in turn were compared with 6,648 similarly genotyped controls recruited at CHOP. The 3,411 cases genotyped on Affymetrix arrays included subjects from the Genetics of Kidneys in Diabetes Study (GoKinD) and the Wellcome Trust Cases Control Consortium (WTCCC) that were then compared with 10,308 similarly genotyped controls, including being derived from non-autoimmune disease related cases from the WTCCC, as well as from the British 1958 Birth Cohort and the UK National Blood Service.

7We compared the power of our meta-analysis to that of the previous largest meta-analysis to date. We have more than double the power of the Barrett et al. meta-analysis to find variants with a relative risk of 1.2 and approximately three times the power to detect variants with a relative risk of 1.1 (Figure S1).

8We used principal components analysis (PCA) in order to minimize the potential impact of population stratification in our case/control sample sets. Eigenstrat 3.0 was employed to remove outliers and to subsequently calculate the principal components in the Illumina and Affymetrix assigned groups separately. The principal components were then used as covariates in a logistic regression, using the software PLINK, to compute the P-values, betas and standard errors which were combined in our fixed effects inverse variance meta-analysis. After controlling for population stratification, the λ in the Affymetrix and Illumina cohorts was 1.11 and 1.17, respectively (see Figure 1 for Q-Q plot). A full description of the correlation of each eigenvector with known continental ancestry appears in Text S1. Mach was used to impute ∼2.54 million SNPs, including HapMap Phase 2 SNPs in the Illumina and Affymetrix datasets in order for the statistics to be uniform and amenable to being combined. Results from the meta-analysis of this resulting ‘discovery’ cohort are shown Table 1 and graphically in Figure 2.

953 SNPs were brought forward to the replication stage based on satisfying the following criteria; however one of these SNPs consistently failed genotyping in the replication effort. The most significantly associated SNP at a given locus if the meta-analysis P-value was less than 1×10−5 (an independent locus was defined as a region for a given focal SNP, where we extended the region in both directions until either 250 kb had been traversed or until reaching another SNP with P<10−5), the Cochran's Q statistic P-value was greater than 0.05 and the locus had not been already reported from a previous GWAS of T1D. A table outlining the results for all previously described T1D associated SNPs plus our strongest associations for known regions associated with the disease are shown in Table 2 and Table S1, respectively. The replication cohort consisted of additional T1D affected trios from the T1DGC and McGill which had not been part of the original discovery cohort. The replication cohort was genotyped using the Sequenom iPLEX system and the results were analyzed using the transmission disequilibrium test in PLINK. Results for both the discovery and replication cohorts for the six SNPs that replicated with P≤0.05 are shown in Table 1 (the full outcomes for all SNPs tested are in Table S2).

10We combined the ‘discovery’ and ‘replication’ meta-analysis P-values using Fisher's combined P-value method implemented in Haploview, comparable to what has been previously employed by others. Three of the SNPs, namely rs539514, rs478222 and rs924043, had combined P-values <5×10−8, the statistical threshold for genome wide significance, while the remaining three, namely rs550448, rs12679857 and rs6547853, failed to reach this bar but were suggestive of association as the alleles yielded both a consistent direction of effect and P-values <0.05 in the replication cohort. These two categories of outcome are summarized in Table 1; in addition, these six SNPs were further investigated with respect to adjustments of the discovery and met-analysis P-values based on the lambdas of each respective cohort (Table S3).

11We have carried out the largest meta-analysis of genome wide genotyped datasets for T1D to date. The replication of three loci using the stratification-free TDT with minimal Mendelian error clearly indicates that they are not false positives due to artifacts such as uncorrected systematic error from stratification or genotyping bias.

12The most significantly associated SNP (rs539514, P = 5.66×10−11) resides in an intronic region of the LMO7 (LIM domain only 7) gene on 13q22. We investigated the associated region using LocusZoom and determined that it is the only gene residing within the block of linkage disequilibrium harboring the signal (Figure S3). Regional plots showing P-values, linkage disequilibrium, and recombination rate for all SNPs in Table 1 are outlined in the Figures S2, S3, S4, S5, S6, S7. LMO7 encodes a protein that contains multiple domains, including a calponin homology domain, a PDZ domain and a LIM domain. There are multiple LMO7 isoforms already known but their full nature and the actual extent of different isoforms remains unclear. Mice with homozygous deletions of LMO7 display retinal, muscular, and growth retardation. Although the function of LMO7 doesn't clearly relate to the etiology of T1D, LMO7 is expressed in pancreatic islets and thus is a possible biological candidate at this locus; however it should be noted that the retinal, muscular development and islet patterns are a key element in Emery-Dreifuss Muscular Dystrophy, caused by mutations in LMO7, but bears very little similarity to T1D.

13The second most significantly associated SNP (rs478222, P = 3.50×10−9) resides in an intronic region of the EFR3B (protein EFR3 homolog B) gene on 2p23; however the region of linkage disequilibrium is approximately 800 kb and harbors additional multiple genes, including 3NCOA1, C2orf79, CENPO, ADCY3, DNAJC27, POMC, and DNMT3A. (Figure S2). A previous meta-analysis of a subset of the data used in this current study found suggestive association with T1D in the same LD block with the independent SNP, rs2165738(r2 = 0.115), but did not achieve genome wide significance at that time (P = 3.65×10−6); however, we only found modest evidence of association with rs2165738 (P = 4.78×10−3) in our discovery cohort. There has also been association to inflammatory bowel disease height, and BMI reported at this locus, where in both cases the risk allele for increased height or BMI was protective for T1D risk.

14The third most significantly associated SNP (rs924043, P = 8.06×10−9) lies in an intergenic region on 6q27, where the region of association is approximately 900 kb and harbors multiple genes including WDR27, C6orf120, PHF10, TCTE3, C6orf208, LOC154449, DLL1, FAM120B, PSMB1, TBP and PCD2 (Figure S5). In addition, despite not reaching the bar for genome wide significance, we did observe evidence for association at three additional loci (Table 1) containing the candidate genes LOC100128081, TNFRSF11B and FOSL2. Of these, it is notable that TNFRSF11B is a strongly associated locus with bone mineral density, also as a consequence of GWAS,. In addition, the locus harboring LOC100128081 has also been reported in the context of a GWAS of SLE. Further work will be required to fully validate the role of these particular loci in the pathogenesis of T1D.

15The Barrett et al. meta-analysis was able to use British controls with British cases and American controls with American cases. We did not have the same control data to be able to make the same comparisons. In the case of the Affymetrix analysis, some American cases were analyzed with purely British controls and, in the case of the Illumina analysis, some British cases with purely American controls. As such, we were forced to make our corrections using eigenvectors as covariates in our analysis; this will have the effect of modestly weakening the level of significance for associations that vary in allele frequency between the cases and controls, as now the case and controls will both vary with the eigenvectors to some degree. This in effect will make our analysis overly conservative with estimating the true effect of a SNP, and in fact every SNP that had a P-value less than 0.05 in the replication set did indeed have a greater effect than that which was estimated from the discovery set.

16In summary, we provide convincing evidence for the existence of three additional loci associated with the T1D, adding to the repertoire of over 50 loci already demonstrated to be associated with the disease.

17The study was approved by the institutional review board and the ethics committee of each institution. Written informed consent was obtained from each participant in accordance with institutional requirements and the Declaration of Helsinki Principles.

18Cases in the discovery set were obtained from four publically available resources and combined with those from a previous publication for the meta-analysis. Samples descriptions are available on dbgap (http://www.ncbi.nlm.nih.gov/sites/entrez?db=gap) for the T1DGC (http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000180.v1.p1), GoKinD (http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000088.v1.p1), and DCCT-EDIC (http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000086.v2.p1) patients. The WTCCC sample information is available from. Samples from the T1D segment of the WTCCC were used as cases, while controls were derived from the 1958 Birth Cohort, UK Blood Service, Bipolar disorder, Coronary heart disease, Hypertension, and Type 2 Diabetes segments. The remaining cases used in the meta-analysis were previously described.

19The total number of individuals used in the meta-analysis discovery set was 26,890 (9,934 cases/16,956 controls). The replication set consisted of 1120 case-parent trios from the T1DGC and those identified through pediatric diabetes clinics in Canada. The replication set was identical to that used in Hakonarson et al. with an extension of patients identified through pediatric diabetes clinics in Montreal, Toronto, Ottawa, and Winnipeg. All individuals were of Caucasian ancestry. A breakdown of the number of samples in each cohort in the discovery phase and a comparison with the numbers used in the Barrett et al. meta-analysis are shown in Table 3. The minor variation in the number of cases reflects that, despite using slight differences in both quality control and methods for dealing with population stratification, we have comparable numbers of cases from this cohort remaining in our analysis. Primarily, this small difference is due to the fact that we strictly accounted for relatedness and duplicates within and across cohorts in this current setting.

20Power analysis was performed with the genetic analysis calculator which can be found at (http://pngu.mgh.harvard.edu/~purcell/gpc/). Various assumption were made included perfect LD between the causative variant and the markers that were genotyped, an additive genetic model, a disease prevalence of 0.0033 and an alpha of 1×10−5.

21Discovery samples from Philadelphia, Canada, T1DGC, and DCCT-EDIC were genotyped on a mixture of the Illumina HumanHap 550v1, 2, and 3, whereas samples from GoKinD and WTCCC were genotyped on the Affymetrix 500 K Chip. Sequenom iPlex was used to replicate the findings of the meta-analysis in 1,120 affected offspring trios from the T1DGC and from Canada.

22All individuals needed an individual genotyping call rate greater than 0.98 to be included in the analysis pre-imputation and individuals were removed that showed evidence of cryptic relatedness and duplication within and across cohorts using identity-by-state. SNP quality control was performed on all samples pre-imputation. SNPs were excluded from the analysis if the minor allele was below 1%, the genotyping call rate was less than 95%, or the Hardy Weinberg equilibrium P-value was less than 0.00001.

23To control for population stratification, Eigenstrat 3.0 was used to compute the top 10 principal components of the individuals genotyped on the Illumina SNP chips and the Affymetrix SNP chips separately. Individuals were removed from the analysis if they were 6 standard deviations away from the mean of one of the top 10 principal components. After controlling for population stratification, the estimated lambda in the Affymetrix data was 1.11 and 1.17 in the Illumina data.

24Mach 1.0 was used to impute ∼2.54 millions SNPs from the HapMap CEU panel for all individuals. SNPs were excluded after imputation if they had a minor allele frequency less than 0.01 and an r2 value less than 0.3.

25PLINK was used to perform a logistic regression using the 10 principal components as covariates, T1D status as the outcome, and in the case of the Affymetrix cohort, an extra dummy covariate specifying WTCCC or GoKinD cohort membership. Results from the logistic regression of 2,436,110 SNPs from the Affymetrix samples and 2,062,307 SNPs from the Illumina samples separately were combined using inverse-variance meta-analysis in PLINK. A fixed effects meta-analysis was performed and 53 SNPs were chosen for replication who had a fixed effects P-value <0.00001, a Cochran's Q statistic P-value greater than 0.05 and were not previously known to be associated with type 1 diabetes. However one of the SNPs consistently failed during the replication effort.