Predicting Islet Cell Autoimmunity and Type 1 Diabetes: An 8-Year TEDDY Study Progress Report

The Environmental Determinants of Diabetes in the Young (TEDDY) study is a large (N = 8,676) prospective cohort study designed to identify environmental factors influencing or protecting against development of islet autoimmunity (IA) and onset of type 1 diabetes (T1D) (1). This ongoing international study began in 2004. Published work to date has examined IA/T1D risk factors including HLA DR-DQ genotypes, sex, T1D-related single nucleotide polymorphisms (SNPs) (2,3), dietary factors (soluble fiber, probiotics, and infant formula type) (4–6), growth in early life (7), age at autoantibody seroconversion, and the first appearing autoantibody (3). TEDDY-identified risk factors are those factors with statistical significance in time-to-event analyses. However, statistically significant risk factors are not automatically good predictors of disease (8). For example, a risk factor may be associated with or cause a disease but be a poor predictor because it is only one of many causal mechanisms. Conversely, a risk factor may predict a disease and yet not be a component of its pathogenesis. Therefore, quantifying (9) the individual and cumulative/incremental predictive utility of the identified risk factors for IA and T1D is an important part of the investigative focus of the TEDDY study, as it should be in other similarly constructed prospective cohort studies.

The aim of this work is to evaluate the individual and combined ability of TEDDY-identified genetic and environmental risk factors to predict IA and progression to T1D. Specificity, sensitivity and receiver operating characteristic (ROC) curves are fundamental tools for diagnostic test evaluation. Since TEDDY is a prospective cohort study, the status of disease outcome (development of IA or onset of T1D) changes over time and the analysis of the cohort at any given point in time is subject to censoring. Hence, time dependency needs to be incorporated in an assessment of the predictive assessment of these risk factors, which may result in different values of sensitivity and specificity over the course of the study (10,11). In this work, time-dependent analysis of sensitivity, specificity, and ROC curve was constructed to evaluate the identified risk factor for predicting disease outcome at several time points of interest.

TEDDY includes six clinical research centers: Colorado, Georgia/Florida, Washington, Finland, Germany, and Sweden. Written informed consent was obtained for all study participants from a parent or primary caretaker, separately, for genetic screening and participation in the prospective follow-up, beginning at birth. The high-risk genotypes selected for inclusion for participants screened from the general population (GP) were as follows: DRB1*04-DQA1*03-DQB1*03:02/DRB1*03-DQA1*05-DQB1*02:01 (DR3/4), DRB1*04-DQA1*03-DQB1*03:02/DRB1*04-DQA1*03-DQB1*03:02 (DR4/4), DRB1*04-DQA1*03-DQB1*03:02/DRB1*08-DQA1*04-DQB1*04:02 (DR4/8), and DRB1*03-DQA1*05-DQB1*02:01/DRB1*03-DQA1*05-DQB1*02:01 (DR3/3). Additional study inclusion genotypes for first-degree relatives (FDRs) of a subject with T1D were DRB1*04-DQA1*03-DQB1*03:02/DRB1*04-DQA1*03-DQB1*02:02 (DR4/4b), DRB1*04-DQA1*03-DQB1*03:02/DRB1*01- DQA1*01-DQB1*05:01 (DR4/1), DRB1*04-DQA1*03-DQB1*03:02/DRB1*13-DQA1*01-DQB1*06:04 (DR4/13), DRB1*04-DQA1*03-DQB1*03:02/DRB1*09-DQA1*03-DQB1*03:03 (DR4/9), and DRB1*03-DQA1*05-DQB1*02:01/DRB1*09-DQA1*03-DQB1*03:03 (DR3/9). The HLA DR-DQ genotype abbreviations shown in parentheses will be used throughout this article. Genotyping was confirmed by reverse blot hybridization at the central HLA Reference Laboratory at Roche Molecular Systems, Oakland, CA (12), along with the INS-23 Hph1 (rs689), CTLA4 T17A (rs231775), and PTPN22 R620 W (rs2476601) SNP primer pairs. Children enrolled are followed prospectively from 3 months to 15 years, with study visits every 3 months until 4 years of age and every 3 or 6 months thereafter depending on autoantibody positivity. All children who are persistent positive for any autoantibody are followed every 3 months until the age of 15 years or onset of T1D. T1D was defined according to the American Diabetes Association’s criteria for diagnosis. Detailed study design and methods have previously been published (1,13). The study was approved by local institutional review or ethics boards and is monitored by an external evaluation committee formed by the National Institutes of Health.

SNP analysis was performed by the Center for Public Health Genomics at University of Virginia, using the Illumina Immunochip, which is a custom array for genotyping of SNPs selected from regions of the human genome firmly associated with autoimmune diseases (14). The final selection of SNPs containing ∼186,000 SNPs in 186 regions, for 12 autoimmune diseases, was decided by the Immunochip Consortium. TEDDY previously examined whether any of 41 non-HLA SNPs previously shown to be associated with T1D conferred risk for IA (2).

Autoantibodies to insulin (IAA), GAD (GADA), or insulinoma antigen-2 (IA-2A) were measured in two laboratories by radiobinding assays (1,13). In the U.S., all sera were assayed at the Barbara Davis Center for Childhood Diabetes at the University of Colorado Denver; in Europe, all sera were assayed at the University of Bristol, U.K. Both laboratories demonstrated high sensitivity and specificity as well as concordance (15). All positive islet autoantibodies and 5% of negative samples were retested in the other reference laboratory and deemed confirmed if concordant. Persistent IA was defined as confirmed positive autoantibodies to IAA, GADA, or IA-2A in at least two consecutive samples.

For the development of IA, the TEDDY study group has examined HLA DR-DQ genotypes, sex, family history of T1D, T1D-associated non-HLA SNPs presented by the Type 1 Diabetes Genetics Consortium (T1DGC) meta-analysis (2,14), complement genes (16), dietary factors (soluble fiber, probiotic, and infant formula type) (4–6), and growth in early life (6,7).

TEDDY has previously published that diabetes-related persistent confirmed islet autoantibodies first appear singly, and autoantibodies to IAA appear at an earlier age than those to GADA (17). The order of appearance was shown to be related to HLA DR-DQ genotypes. Therefore, TEDDY also examined all of the aforementioned IA-associated risk factors, along with early infection-related conditions, on the risks of developing IAA or GADA as the first appearing indication of autoimmunity (IAA appearing first [IAA-first] or GADA-first, respectively) using multivariate cause-specific proportional hazards modeling (18).

For T1D, the analysis included age at multiple autoantibodies, more than a single autoantibody as the first appearing indication of seroconversion, and female sex along with SNPs rs10517086, rs1534422, and rs2327832 in TNFAIP3, which were associated with an increased risk of progression to T1D, and SNP rs1004446 in INS, which was associated with decreased risk (3).

For each outcome (IA, IAA-first, GADA-first, or progression to T1D), we evaluated individual predictive ability of the TEDDY-identified risk factors by calculating time-dependent sensitivity (Se) and specificity (Sp) at 6 years of age (see Supplementary Data for 2 and 10 years follow-up). Family history of T1D was categorized as FDR versus GP, HLA genotypes were categorized as DR3/4 versus others, type of formula was categorized as extensively hydrolyzed cow’s milk–based formula versus others, each of the SNPs was categorized as major (no copy of minor allele) versus minor/heterogeneous (one or two copies of minor alleles), and first appearing indication of seroconversion was categorized as more than a single autoantibody versus others. Continuous risk factors, e.g., the weight z score at 1 year and age at multiple autoantibodies onset, were evaluated using median value as the cutoff.

Since we were also interested in the incremental predictive ability of each risk factor, we performed a forward selection procedure, with an entry level (P < 0.1) to order and include risk factors from all of the identified risk factors by Cox proportional hazards regression. Then, we fitted Cox models by adding one risk factor at a time according to the order determined by the forward selection. For each of the fitted Cox models, a risk score was defined as the estimated linear predictor in the log hazard function. A time-dependent ROC curve at each of the specified time points was calculated for each of the risk scores (prediction models). The Youden’s index (J), Se + Sp − 1 (9), and area under the curve of the time-dependent ROC (AUC) were calculated. For the analyses presented herein, J is defined at its maximum. The AUC quantifies the probability that the risk scores from a randomly selected pair of diseased and nondiseased individuals are correctly ordered.

The forward selection procedure was performed using the SAS software (version 9.4; SAS Institute, Cary, NC). The package “survival” (19) was used to fit Cox models and the package “timeROC” (20) was used to calculate time-dependent sensitivity, specificity, and ROC in R 3.4.2 (R Foundation for Statistical Computing, Vienna, Austria).

As of 31 July 2017, 7,777 of 8,676 (89.6%) enrolled children who were HLA DR-DQ genotype eligible and had Immunochip data available were included in the analyses. The median age of follow-up for the 7,777 children was 9.1 years (interquartile range 7.9–10.4). As summarized in Table 1, a total of 736 (9.5%) children developed one or more persistent autoantibodies (GADA, IAA, or IA-2A), and among them 281 (38%) children had IAA only as the first appearing autoantibody (IAA-first), 316 (43%) children had GADA only as the first appearing autoantibody (GADA-first), 17 (2%) children had IA-2A only as the first appearing autoantibody (IA-2A-first), and 104 (14%) children had two autoantibodies and 18 (2%) children had three autoantibodies appearing simultaneously. Altogether, 434 subjects developed two or more autoantibodies, and of these, 219 (50.5%) children developed T1D.

The sensitivity (Se), specificity (Sp), and Youden’s index (J) of each risk factor for predicting IA, IAA-first, and GADA-first at 6 years are presented in Table 2. Among all the identified risk factors, HLA genotype (DR3/4 vs. others) is the best predictor for IA (J = 0.117) and SNP rs2476601 (PTPN22) is the best predictor for IAA-first (J = 0.123), whereas family history of T1D (J = 0.100 and 0.111) for both, SNP rs2476601 for IA (J = 0.112), and HLA genotype for IAA-first (J = 0.104) achieved at least 0.10 in Youden’s index. For GADA-first, weight at 1 year is the best predictor (J = 0.114), followed by HLA genotype (J = 0.103) and SNP rs2292239 in ERBB3 (J = 0.112). SNP rs3184504 in SH2B3 also achieved >0.10 in Youden’s index (J = 0.105). Similar results have been observed in predicting IA and IAA-first, while SNP rs2476601 is the best predictor for GADA-first (J = 0.191) at 2 years (Supplementary Table 1). For prediction of IA, IAA-first, and GADA-first at 10 years (Supplementary Table 2), HLA genotype is the best predictor, with Youden’s index of J = 0.105, 0.113, and 0.09, respectively. No other risk factor achieved >0.10 in Youden’s index. From the plots of (Se, 1 − Sp) in Fig. 1, it is easy to see that these risk factors individually lie close to the diagonal line, which is interpreted as providing little predictive value for IA (IA, IAA-first, or GADA-first) or progression to T1D.

In the forward selection for IA, a total of 12 risk factors were added into the prediction models in the order of family history of T1D, HLA DR-DQ genotype, rs2476601 (PTPN22), rs3184504 (SH2B3), rs1004446 (INS), weight z score at 1 year of age, rs2292239 (ERBB3), rs1270876 (CLEC16A), sex, rs4948088 (COBL), probiotics exposure before 28 days, and type of first formula introduced before 3 months of age. The AUC and the Youden’s index (J) of each prediction model for predicting IA at 6 years are given in Table 3. The AUC was 0.552 (95% CI 0.534, 0.569) with the inclusion of family history of T1D and increased to 0.678 (95% CI 0.655, 0.701) with the full model at 6 years. Both AUC and Youden’s index increased when additional risk factors were added into the prediction models. The incremental AUC (∆AUC) and Youden’s index (∆J) introduced by each newly added risk factor are also given in Table 3. The first four risk factors (family history of T1D, HLA genotype, SNP rs2476601 in PTPN22, and SNP rs3184504 in SH2B3) contributed the most. Similar results have been observed in predicting IA at 2 and 10 years (Supplementary Tables 3 and 4 and Fig. 1A). The ROC curves to predict IA at 2, 6, and 10 years using the prediction model with all of the 12 risk factors are plotted in Fig. 2A. As we can see, the ability to predict IA gets worse as the time frame gets lengthened.

For IAA-first, the AUCs at years 2, 6, and 10 reached 0.704 (95% CI 0.665, 0.744), 0.707 (95% CI 0.676, 0.739), and 0.683 (95% CI 0.647, 0.719), respectively, when all eight identified risk factors were combined into the prediction model (Table 3, Supplementary Tables 3 and 4, and Fig. 2B). The first five added risk factors (family history of T1D, HLA DR-DQ genotype, SNP rs2476601 in PTPN22, SNP rs1004446 in INS, and SNP rs3184504 in SH2B3) contributed the most in terms of the AUC. For GADA-first, the AUCs at years 2, 6, and 10 reached 0.683 (95% CI 0.618, 0.747), 0.686 (95% CI 0.651, 0.722), and 0.639 (95% CI 0.603, 0.674), respectively, when all risk factors were combined into the prediction model (Table 3, Fig. 2C, and Supplementary Tables 3 and 4).

The sensitivity and specificity of each risk factor for predicting progression to T1D at 3 years after onset of multiple autoantibodies are presented in Table 2 and are plotted in Fig. 1D. Among all identified risk factors, age at onset of multiple autoantibodies is the best predictor for progression to T1D at 3 years in terms of Youden’s index (J = 0.228). First appearing indication of seroconversion ranked second in terms of Youden’s index (J = 0.127). Similar results have been observed in predicting T1D at 1 and 5 years (Supplementary Tables 1 and 2 and Supplementary Figs. 1D and 2D).

In the forward selection, all seven risk factors were added into the predictive model in the order of age at onset of multiple antibodies, first appearing indication of seroconversion, SNP rs1534422, sex, and SNPs rs2327832, rs1004446, and rs10517086. The AUC and Youden’s index of each prediction model for predicting progression to T1D at 3 years are given in Table 3 and at 1 and 5 years are given in Supplementary Tables 1 and 2. Age at onset of multiple antibodies is a strong predictor for progression to T1D, and other risk factors contributed only a little in additional predictive value. The AUCs of the prediction model combining all identified risk factors reached 0.781 (95% CI 0.716, 0.845), 0.706 (95% CI 0.649, 0.762), and 0.660 (95% CI 0.600, 0.720), respectively, to predict progression to T1D at years 1, 3, and 5 (Fig. 2D). As to the ability to predict IA, the prediction models gave a better prediction of progression to T1D at an earlier time point compared with a later time point.

While TEDDY has identified genetic and environmental markers significantly associated with the risk of IA and T1D, this work has evaluated their individual and combined contribution to predicting IA and progression to T1D based on time-dependent specificity, sensitivity, and the corresponding ROC curve. Despite the strong statistical associations, the results described herein show that these associations are relatively weak predictors. The Youden indices are generally <0.4 and the AUC <0.8, which put the predictors no better than fair. Interestingly, the ability to predict IA or T1D decreases with an increasing time horizon. In addition, the ability to predict IAA-first and GADA-first is no better than IA, suggesting the strong influence of age (most IAA-first occurs at an early age compared with GADA-first, which is consistent with our previous results [17,18]).

With respect to prediction of IAA-first, HLA genotype contributed the most toward prediction after family history of T1D. The contribution of three other genetic factors (rs2476601 in PTPN22, rs1004446 in INS, and rs3184504 in SH2B3) contributed equally, raising the J index from 0.115 to 0.316 with the others contributing little improvement. With respect to predicting GADA-first, rs3184504 in SH2B3_T contributed the most followed by rs2476601 in PTPN22 and family history of T1D contributing equally. After consideration of these risk factors, additional contribution to prediction from the others was really nominal. Age at onset of multiple autoantibodies was the most significant predictor of T1D, with the other risk factors also making a very nominal contribution. Of note, the TEDDY cohort is young, and these results may be affected by the early seroconversion to IAA-first autoantibody positivity, which diminishes as age increases. Additional follow-up may increase the ability to predict T1D following GADA-first seroconversion.

The ability to estimate sensitivity and specificity, in a time-until-event setting, provides a means by which to evaluate the contribution of statistically significant risk factors identified by using Cox proportional hazards models. Those models describe contributions to risk estimates (hazard rates) that are only the positive predictive value. The methods included herein provide an ability to assess the extent to which risk factors can discriminate between those who will and will not get disease.

The incorporation of significant risk factors identified from TEDDY’s nested case-control study, such as the observation of lower IA odds ratio associated with higher average 25(OH)D concentration over all time points using conditional logistic regression analysis (21), into the prediction modeling requires analytic methods that extend the results to the entire cohort. The matched concordance index (C-index) (22) of 0.56 (95% CI 0.53–0.59) suggests that it adds little to the data presented herein.

The factors identified and the analysis presented in this article are not without their limitations. The TEDDY cohort is relatively young, with a selected high-risk HLA from selected geographic settings. The results reported herein may not remain the same when the cohort ages or if applied to other cohorts. But this point speaks to generalizability of each of the risk factors, not to the methods described here, which provide a means to assess their relative contribution to prediction in the same cohort from which they were identified. It is important to have a validation data set, and we accept this limitation in TEDDY overall. Finally, the analytic methodology that we used to combine the risk factors has its own limitations, yet follows accepted practice for such an analysis (23).

As additional risk factors are identified and prediction improves, the results can provide insights into causal mechanisms. TEDDY is uniquely poised to make this contribution as it adds environmental exposures, pathogens, and extensive genomics to the risk factors identified thus far. Going forward, the presented evaluation framework will serve as a summary tool to evaluate these new risk factors as they are identified in TEDDY.

Acknowledgments. The authors thank Sarah Austin-Gonzalez with the Health Informatics Institute at the University of South Florida for assistance with preparing the figures.

Funding. The TEDDY Study is funded by the National Institute of Diabetes and Digestive and Kidney Diseases, National Institute of Allergy and Infectious Diseases, National Institute of Child Health and Human Development, National Institute of Environmental Health Sciences, JDRF, and the Centers for Disease Control and Prevention (grants U01 DK63829, U01 DK63861, U01 DK63821, U01 DK63865, U01 DK63863, U01 DK63836, U01 DK63790, UC4 DK63829, UC4 DK63861, UC4 DK63821, UC4 DK63865, UC4 DK63863, UC4 DK63836, UC4 DK95300, UC4 DK100238, UC4 DK106955, UC4 DK112243, and UC4 DK117483, and contract no. HHSN267200700014C). This work was supported in part by National Institutes of Health/National Center for Advancing Translational Sciences Clinical and Translational Science Awards to the University of Florida (UL1 TR000064) and the University of Colorado (UL1 TR001082).

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. J.P.K., X.L., K.V., B.A., W.A.H., M.J.R., J.-X.S., J.T., A.-G.Z., and Å.L. attest to meeting ICMJE (International Committee of Medical Journal Editors) uniform requirements for authorship by making substantial contributions to conception and design of this manuscript; acquisition, analysis, and interpretation of the data; drafting or revising the manuscript for intellectual content; and giving final approval of the published version. J.P.K. designed the study, proposed the analysis, interpreted the findings, and wrote the manuscript. X.L. performed the analysis and drafted/revised the manuscript. K.V. provided input on the analytical plan and interpretation of the results and drafted/revised the manuscript. B.A., W.A.H., M.J.R., J.-X.S., J.T., A.-G.Z., and Å.L. designed the study and reviewed/edited the manuscript. J.P.K. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Prior Presentation. Parts of this study were presented in abstract form at the 78th Scientific Sessions of the American Diabetes Association, Orlando, FL, 22–26 June 2018.