Risk of Type 1 Diabetes Progression in Islet Autoantibody-Positive Children Can Be Further Stratified Using Expression Patterns of Multiple Genes Implicated in Peripheral Blood Lymphocyte Activation and Function

Type 1 diabetes (T1D) is a chronic autoimmune disease resulting from the targeted destruction of insulin-secreting pancreatic islet β-cells. Islet autoantibodies, markers of active islet autoimmunity, can be detected years, and even decades, before the appearance of clinical symptoms (1). The long asymptomatic period between the appearance of islet autoantibodies and disease onset provides a window of opportunity for T1D prevention in subjects who are autoantibody positive (AbP). At least one of four major islet autoantibodies (insulin autoantibody, GAD antibody, IA-2 antibody, and zinc transporter 8 antibody) is detected in >90% of patients with newly diagnosed T1D (2). These autoantibodies have become the gold standard for identifying at-risk subjects from first-degree relatives (FDRs) of T1D patients as well as the general population (3). Further improvement of risk prediction using autoantibodies has clear scientific and clinical value. Subjects with multiple islet autoantibodies have a high projected risk within 10 years (69.7% [95% CI, 65.1–74.3%]), whereas the presence of a single autoantibody shows a low risk (14.5% [95% CI, 10.3–18.7%]) (4). Furthermore, AbP subjects have variable progression to T1D, with a prediabetes period ranging from 0 to 20 years. Given the variable length of time, further stratifying these individuals for more accurate prediction to clinical disease would be advantageous. Although age at seroconversion and titer of autoantibodies can further improve risk prediction (4,5), additional biomarkers are still needed. Considerable efforts have been devoted to the development of genetic and metabolic biomarkers based on AbP prospective cohorts or T1D prevention trials. HLA (4,6–8), non-HLA (9–12) genetic markers, and metabolic risk scores (13–15) have shown certain levels of improvement for stratifying the risk of AbP subjects. However, the practical potential of these markers is limited by either their intrinsic deficiencies or low predictive values.

Gene expression profiles are expected to dynamically change during disease progression and treatment. Therefore, gene expression patterns may serve as potential biomarkers for risk stratification and therapeutic monitoring. Several studies have examined gene expression changes related to T1D and identified a large number of genes that may differ in expression levels among healthy control subjects, AbP subjects, and T1D patients (16–22). However, these studies have been limited by their cross-sectional design and, hence, hardly suggest biomarker potential. The present study identified five genes that in combination can serve as biomarkers to stratify progression risk in AbP subjects. Our strategy first used microarray data to discover gene expression changes associated with differential progression from AbP to T1D and then validated the top-27 genes using quantitative RT-PCR data from independent AbP subjects from the Diabetes Autoimmunity Study in the Young (DAISY) cohort.

A total of 104 AbP subjects who were consecutively seen in DAISY and followed until February 2012 were included in the analyses. AbP status was identified based on the presence of at least one of the following three autoantibodies: insulin autoantibody, GAD antibody, and IA-2 antibody. By the cutoff date, diabetes developed in 39 of the 104 AbP subjects, with a median follow-up time (from first AbP) of 5.64 years. The median follow-up time of the 65 nonprogressors was 8.9 years. Diabetes was diagnosed according to American Diabetes Association criteria. Demographic information on the distribution of age and sex, age at the appearance of first antibody, FDR status, genetic risk [classified by HLA genotype (23)], number of antibodies, and follow-up time (after first AbP) is summarized in Table 1. The 104 subjects were split into two study phases: discovery (microarray) and validation (real-time RT-PCR). Thirty-six subjects were selected for the discovery phase, with progressors (n = 21) and nonprogressors (n = 15) matched for age, sex, age at first AbP, FDR status, genetic risk, and number of autoantibodies (most with two or more AbP). The remaining 68 subjects (18 progressors and 50 nonprogressors) were included in the validation phase. Compared with nonprogressors, the progressors in the validation data set were younger, had an earlier age of first antibody appearance, had a higher frequency of multiple autoantibodies, and had more frequently FDRs with diabetes. Demographic information for the subjects in the two data sets is summarized in Table 1. A dot plot of the distribution of age of first AbP, age of sample collection, and age of T1D/last visit for each subject is shown in Supplementary Fig. 1. For most subjects, gene expression was analyzed at a relatively early time point following first antibody detection, with a mean interval time of 2.01 ± 2.17 years.

Peripheral blood (2.5 mL) was immediately preserved in PAXgene RNA tubes. After 2 h at room temperature, the tubes were stored at −80°C, and total RNA was extracted within a few weeks using a QIAGEN kit specially developed for the PAXgene RNA tubes. RNA quality and quantity was determined using NanoDrop ND-1000. Extracted RNA was stored at −80°C until use.

Gene expression profiling was performed using the HumanRef-8 Expression BeadChip (Illumina, San Diego, CA). The BeadChip targets ∼24,526 probes. Probe preparation and hybridization were carried out according to manufacturer’s recommendations. The microarray data are Minimum Information About a Microarray Experiment compliant and have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus repository and are accessible through Gene Expression Omnibus Series accession number GSE51058.

High-throughput real-time RT-PCR was performed on the BioMark 96.96 Dynamic Array (Fluidigm) with TaqMan Gene Expression Assays (Applied Biosystems). Briefly, an aliquot of total RNA was converted to cDNA through a High Capacity cDNA Reverse Transcription Kit (Applied Biosystems), then cDNA was preamplified with pooled gene assays. Preamplified cDNA was loaded into BioMark Dynamic Array chips using the NanoFlex IFC controller. Threshold cycle, as a measurement of relative fluorescence intensity, was extracted by BioMark Real-Time PCR Analysis software. Threshold cycle values were directly used in data analysis after normalization to eight preselected reference genes: RNF31, TM9SF4, ABL, PPIB, ESD, FPGS, MRPL19, and TRAP1 (Supplementary Table 1). Technical reproducibility was examined using a pooled cDNA sample aliquot.

All analyses were performed using the R statistics package and environment for statistical computing (www.r-project.org). The lumi package in R was used for preprocessing and normalization of the microarray data using variance stabilizing transformation and robust spline normalization. Probes not detected in >50% of the samples were excluded from further analysis. We used Cox proportional hazard models to evaluate the impact of gene expression levels on diabetes-free survival. Diabetes-free survival was calculated as time from first antibody-positive date to the T1D diagnosis date. Subjects with no evidence of disease were censored at the date of the last follow-up visit. The subjects were divided into two groups (low and high) using the median gene expression value of each individual gene as a cutoff. Kaplan-Meier plots were used to compare differences in diabetes-free survival between these groups. A combination of log-rank test P value and absolute value of the hazard ratio (HR) was used for selecting significant genes. Only genes having a significant effect in the univariate analysis were entered into the multivariate analyses to assess the combined effect of various genes on diabetes-free survival. For the models using multiple genes, self-organizing map cluster analysis (24) was used to divide the subjects into two clusters (C1 and C2) based on the gene expression patterns. Kaplan-Meier plots were used to compare the differences in diabetes-free survival between these groups. The consistency between models was assessed by calculating the percentage of subjects classified in the same group by model. Cox regression analysis was used to evaluate the independent prognostic value of gene expression by including known confounding factors, such as age at first AbP, the number of autoantibodies present, genetic risk, and FDR status, in the regression models. The Akaike information criterion (AIC) was calculated to compare the quality of fit for each model.

Microarray analysis was carried out on peripheral blood mononuclear cells from 36 AbP subjects, 21 of whom had progressed to T1D by the end of follow-up. Among the 12,904 genes tested, 454 showed significant associations (P < 0.05) with T1D progression, and 112 of them had an HR >3 (or <0.33).

To validate these findings, the top-27 genes were selected for real-time RT-PCR confirmation in an independent cohort that included 68 AbP subjects (18 progressors and 50 nonprogressors). Genes were first selected based on HR and P values (HR or 1/HR close to 3, P < 0.02) and were further prioritized based on functional relevance to autoimmunity (especially T1D), inflammation, cellular activation, or proliferation. Three randomly selected negative control genes (XRCC2, FKBP11, and ORM1) were also included for validation (Table 2). Among the 27 candidate genes, 3 (BACH2, TXNDC5, and CDC20) had a confirmed high HR (3.9, 3.3, and 2.8, respectively) and significant P values (P < 0.05), whereas 2 (IGLL3 and EIF3A) had a confirmed good HR (2.6 and 2.3, respectively) and marginal P values (P = 0.055 and 0.074, respectively). The remaining 22 candidate genes and the 3 negative control genes were excluded from further study because of insignificant P values (P > 0.1) or inconsistent HRs in the discovery and validation data sets (Table 2). As shown in Fig. 1, the higher expression of these five genes is associated with a higher risk of T1D progression. The difference in mean gene expression levels of the high-expression subset and low-expression subset are also shown in Fig. 1 (fold change 1.6–3.9 in the validation data set) and are consistent with the microarray data set for four of the five genes.

It is well known that predictive power using multiple genes may be significantly better than that of any single gene. Thus, we evaluated the predictive value of models containing combinations of the five genes selected by the single-gene analysis (BACH2, IGLL3, TXNRD5, CDC20, and EIF3A). With five genes, there is a total of 25 two-, three-, and four-gene models. In the discovery data set, 17 of the 25 models can distinguish two clusters of subjects with significant differences in T1D progression (Supplementary Table 2). Further assessment of these models in the validation cohort confirmed 14 of the 17 models to have predictive value. It should be noted that not every confirmed combination model had an increased predictive power over the single genes. Four models were preferred for their larger HR values, more significant P values, and lower AIC values (Table 3). Model BI (comprising the BACH2 and IGLL3 genes) had an HR of 7.34 (adjusted P = 0.004) and AIC of 122.02. Model ICE (comprising IGLL3, CDC20, and EIF3A) had an HR of 6.32 (adjusted P = 0.007) and AIC of 122.02. Model BICE (comprising BACH2, IGLL3, CDC20, and EIF3A) had an HR of 7.76 (adjusted P = 0.003) and AIC of 121.12. Finally, model BITE (comprising BACH2, IGLL3, TXNDC5, and EIF3A) had an HR of 7.76 (adjusted P = 0.003) and AIC of 121.12. As shown in Fig. 2, these four models consistently stratify AbP subjects with differential progression rates in both the discovery and the validation data sets.

Having shown the utility of the multigene models in risk stratification of AbP subjects, we next assessed the consistency of risk stratification for individual subjects using the various models. Fig. 3 shows the consistency among all possible combinations of these top-four models. A major conclusion from these results is that 75% of the individuals are consistently classified into the same progression group by all four models. Further examination of the data suggests that model ICE has less consistency than the other three models (BI, BICE, and BITE). The mean consistency among models excluding ICE is >89% for both data sets, whereas the mean consistencies for models including ICE are 82.6% and 84.4% for the discovery and validation data sets, respectively. Indeed, ICE has the lowest HR and highest AIC in both the discovery and the validation data sets compared with the other three models (Fig. 2, Table 3).

Several demographic parameters and molecules are well-established risk factors associated with T1D progression. These factors include age at first appearance of autoantibodies, the number of autoantibodies present, genetic factors, and FDR status. With the exception of FDR status, significant differences in these factors were observed between the progressors and the nonprogressors (Table 1). To examine whether these known factors were confounding factors for the gene expression models on T1D risk stratification, Cox regression analysis was used to analyze both the discovery and the validation data sets. According to multivariate analyses, three of the four gene expression models (BI, BICE, and BITE) are independent predictors of T1D progression with significant adjusted P values (P < 0.05) (Table 4).

The study also suggests that using gene expression in conjunction with other known risk markers improves the risk stratification for AbP subjects. Consistent with previous findings, subjects with more than one autoantibody have a significantly higher risk for T1D progression than those with only one autoantibody (HR 3.07, P < 0.05) (Supplementary Fig. 2A). Of note, the subjects with more than one autoantibody can be further stratified into a high-risk and a low-risk group by model BI (HR 7.12, P < 0.05) (Supplementary Fig. 2B).

Islet autoantibodies are the current gold standard for risk stratification for the development of T1D and have been used in large population-based studies, including TEDDY (The Environmental Determinants of Diabetes in the Young) and many clinical trials for T1D prevention (25–27). However, testing for islet autoantibodies is not specific enough even when multiple autoantibodies are detected. The aim of the current study was to investigate whether gene expression patterns can be used as biomarkers to complement islet autoantibodies for T1D risk stratification. We demonstrate that expression changes of certain genes are associated with T1D progression from autoimmunity. The predictive value of single genes is generally moderate (HRs of 2.4–3.9) and may be further limited because the risk assignment by different genes may not be consistent enough. Several multigene models have shown much better predictive power than individual genes. One two-gene model (BACH2 and IGLL3) identified two AbP subgroups with very different progression rates (HR 7.3). Adding two other genes (CDC20 and EIF3A or TXNDC5 and EIF3A) did not significantly improve the HR values (HR 7.8). However, the two four-gene models are highly consistent in risk assignment (consistency of 100% in the validation data set). Because risk of T1D progression among AbP subjects is influenced by many other factors, such as family history of diabetes, genetic factors, and number of autoantibodies, it was necessary to determine whether the risk associated with differential gene expression is independent of these known factors. Cox regression analyses suggested that risk attributed to gene expression could not be accounted for by the risk factors examined in this study. Indeed, gene expression can be used together with the number of autoantibodies to improve risk stratification.

This study has demonstrated for the first time in our knowledge that gene expression patterns may serve as biomarkers to stratify T1D risk among AbP children. When compared with genetic markers (HLA and non-HLA) (4,6–12) and metabolic markers (13–15), the present gene expression markers have several strengths. First, the intrinsic deficiencies of genetic and metabolic markers limit their usefulness for T1D prevention. For example, genetic markers cannot be used as surrogate markers for therapeutic outcomes. Whereas metabolic anomalies are indicative of pancreatic islet cell dysfunction, they are only useful markers at a very late stage of disease progression when the opportunity for disease prevention may have passed. In contrast, gene expression markers can be used for both early disease prediction and therapeutic evaluation because of their dynamic change during disease progression. Second, the identified gene expression markers were shown to further stratify AbP subjects into subgroups with different progression rates to T1D. High HR values (up to 7.76), consistency of their predictive values between discovery and validation data sets, and high model consistency (75–100%) indicate that these gene expression biomarkers have a potentially strong predictive power and reliability. Finally, cohorts in most previous studies were restricted to children with a family history of T1D, such as the genetic studies in BABYDIAB (10–12), the Belgian Diabetes Registry (7), and DiMe (Childhood Diabetes in Finland) (8) and metabolic marker studies in Diabetes Prevention Trial–Type 1 (13–15). The DAISY cohort used in the current study includes children from the general population and with T1D relatives. The Cox regression results suggest that the gene expression markers are independent of family history.

The limitations of this study include the absence of zinc transporter 8 antibody information and the relatively small number of AbP subjects, especially when split into two sample sets. Despite the potential utility of these gene expression models in risk stratification, it should be noted that as a result of population differences and sampling bias, the best genes and models may not be selected in this study. However, the experimental and analytical approaches advocated in this study should serve as a model for future studies using much larger sample sizes. Also essential is that the candidate genes and models identified in this study are further validated in other cohorts.

From a functional perspective, the genes identified in this study may provide novel insights into the disease mechanism. BACH2 and IGLL3 are two B-cell–related genes. IGLL3 encodes the Igλ–like polypeptide 3, which is expressed mainly in pre-B cells and involved in B-cell development (28). The present study is the first in our knowledge to suggest a potential role of this gene in T1D or other autoimmune diseases and should be investigated further in T1D human and animal models. BACH2 is a B-cell–specific transcription factor critical for Ig class switching. B cells with lower levels of BACH2 mature directly to IgM plasma cells, whereas cells with higher levels undergo class switching and differentiate to other cell types, mainly IgG-secreting plasma cells, including autoantibody-secreting cells (29). Higher expression of BACH2 is associated with a faster and higher progression rate from AbP to T1D, consistent with the concept that antibody switching may be a potential risk factor that can differentiate progression risks from AbP to T1D (30). Genetic polymorphisms within the BACH2 gene are associated with several autoimmune diseases, including T1D (31), Crohn disease (32), celiac disease (33), multiple sclerosis (34), and vitiligo (35). Although T1D was considered an ultimately T-cell–mediated autoimmune disorder, accumulating evidence in animal models and human studies indicates a crucial role for B cells in T1D development. The present findings further strengthen the connection between B cells and T1D.

EIF3A and CDC20 are involved in cell cycle regulation (36,37). The association of the two cell cycle–related genes with T1D progression is consistent with a suspected increase in the activation and proliferation of immune cells, especially B cells and T cells, in AbP subjects who are progressing to T1D. TXNDC5 encodes the thioredoxin domain-containing 5 protein, a protein-disulfide isomerase that possesses thioredoxin activity. TXNDC5 may decrease receptor binding activity of insulin by catalyzing the reduction of insulin disulfide bonds (38). TXNDC5 may increase T1D progression risk by promoting glucose intolerance through reducing insulin receptor binding activity, a hypothesis that should be verified in future studies.

In summary, this retrospective study using transcriptomic profiling and RT-PCR validation suggests that five genes may be involved in T1D pathogenesis and may serve as biomarkers for risk stratification among AbP subjects. The findings further indicate that gene expression patterns of multiple genes have much better predictive value and higher consistency in risk assessment than individual genes.

Funding. This work was supported by grants from the National Institutes of Health (R01-DK-32083, DK-32493, R01-DK-050979) to M.R. and (4R33-DK-069878 and 2R01-HD-37800) to J.-X.S. and from the JDRF (17-2013-535) to M.R.

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

Author Contributions. Y.J. contributed to research design, data acquisition, and manuscript preparation. A.S. and S.B. contributed to the statistical analysis and data interpretation. C.D. and H.L. contributed to data generation. D.H. and K.B. contributed to sample and clinical information collection. M.R. contributed to the clinical data and samples, data interpretation, and manuscript preparation. J.-X.S. supervised the overall study design and data interpretation and edited the manuscript. J.-X.S. is the guarantor of this work and, as such, had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.