A Regulatory T-Cell Gene Signature Is a Specific and Sensitive Biomarker to Identify Children With New-Onset Type 1 Diabetes

The precise mechanisms underlying the failure to control autoreactive T-cell expansion in type 1 diabetes (T1D) are unclear, but evidence suggests that changes in FOXP3⁺ regulatory T cells (Tregs) may be involved. People lacking functional Tregs due to genetic mutation are born with T1D (1), and in NOD mice, Treg deletion accelerates T1D progression, whereas transfer of Tregs prevents disease (2). In patients with T1D, Treg proportions in blood appear to be normal (3), although accurate enumeration is difficult because activated conventional T cells (Tconvs) also express FOXP3 and CD25. Some evidence points to impaired Treg function in T1D (4–7), possibly related to differentiation into cells producing interleukin (IL)-17 (6) or interferon-γ (7), but in vitro suppression assays may not reflect in vivo function because cells are removed from their inflammatory environment. T1D Tregs may also have diminished IL-2–mediated activation of STAT5 (8), resulting in destabilized FOXP3 expression and diminished Treg stability.

Although many clinical trials in T1D have aimed to enhance Tregs, the lack of a simple method to accurately monitor Tregs in blood is a major limitation that hinders identification of subjects with dysfunctional Tregs and, thus, the tracking of changes in Tregs throughout disease progression. We report a new gene signature of human Tregs created by comparing the gene expression profiles of memory and naive Tregs with those of memory and naive Tconvs ex vivo and after in vitro activation. Measurement of the Treg signature in T1D revealed a subset of genes that accurately classifies children with new-onset T1D versus healthy children. These data provide new evidence that T1D Tregs have intrinsic changes and provide a clinically applicable way to monitor for changes in Tregs in T1D.

Peripheral blood from 12 adult and 24 pediatric control subjects, 29 patients with new-onset T1D, and 27 patients with established T1D was obtained upon written consent and age-appropriate assent as approved by The University of British Columbia Clinical Research Ethics Board. Samples were collected on average 1.5 (range 0.9–3.4) and 28.9 (9.4–50.2) months after diagnosis of new-onset and established T1D, respectively.

Adult blood was pre-enriched for CD25⁺ and CD25⁻ cells (magnetic bead sorting; Miltenyi) and sorted into naive (CD4⁺CD25^hiCD45RA⁺) and memory (CD4⁺CD25^hiCD45RA⁻) Tregs and naive (CD4⁺CD25⁻CD45RA⁺) and memory (CD4⁺CD25⁻CD45RA⁻) Tconvs. RNA was extracted directly or after 40-h stimulation with α-CD3/α-CD28–coated beads (Life Technologies; 1:1 bead:cell ratio, 100 units/mL IL-2 [Chiron] in X-VIVO 15 containing 5% human serum [both Cambrex] and penicillin/streptomycin and GlutaMAX [both Life Technologies]). To test suppressive capacity, CD4⁺ T cells were stimulated (1:8 bead:cell α-CD3/α-CD28–coated beads) with and without Treg or Tconv subsets for 6 days, and wells were pulsed with [3H]-thymidine (1 μCi/well; Amersham) for the last 16 h.

For some adult and all pediatric cryopreserved peripheral blood mononuclear cells (PBMCs), CD4⁺ T cells were enriched (STEMCELL Technologies), sorted into Tregs (CD4⁺CD25^hiCD127^lo) and Tconvs (CD4⁺CD25^loCD127^hi), and lysed in RLT lysis buffer (QIAGEN). Where possible, aliquots were stained for FOXP3 (allophycocyanin; eBiosciences) to assess purity. For some pediatric samples (n = 6–8), unfractionated PBMCs were lysed in RLT lysis buffer for NanoString.

DNA was extracted from male Tregs and Tconvs, bisulphite converted (EZ DNA Methylation-Direct Kit; Zymo Research), and assayed for Treg-specific demethylated region (TSDR) methylation using the PyroMark Gold Q96 PCR Kit (QIAGEN) (primers: EpigenDx) on a PyroMark Q96 ID pyrosequencer (QIAGEN).

mRNA was extracted, reverse transcribed, and amplified, and cDNA was hybridized to Affymetrix human genome U133A 2.0 microarrays (GEO accession GSE76598). After preprocessing, normalization, and quality control, four samples were excluded. The data were analyzed with the Marker Selection module (GENE-E microarray analysis software) with a false discovery rate cutoff of <25%. Two rounds of Gene Set Enrichment Analysis were performed with and without a cutoff (750-count difference between Tregs and Tconvs) and including publicly available data sets of mouse (9,10) or human (11–15) Treg/Tconv gene expression profiles.

mRNA expression was measured with a custom nCounter reporter probe set (31 genes of interest and 6 reference genes) by using 100 ng total RNA or 2–4 μL whole-cell lysate (5,000 cells/μL lysis buffer). Before biomarker discovery analysis or hierarchical clustering, batch correction (ComBat [16]), positive control sum normalization, background subtraction (with counts <1 changed to 1), sum normalization, and log₂-transformation were performed (Bioconductor [17]).

The hierarchical clustering module (public GenePattern server, Broad Institute [18]) was used to generate heat maps and cluster trees by pairwise average linkage according to Pearson correlation coefficient.

Principal component analysis (PCA) was performed on normalized data in R (19). Using elastic net calculations and 200- × 5-fold cross-validation, the 30-gene panel (1 gene was undetectable and excluded) performance was tested (α = 0) by comparing disease groups per cell subset. Classification performance was tested in R (19) by using area under the receiver operating characteristic curve (AUC), sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and error rate. Biomarker discovery was repeated with α = 0.65 (mean panel size 10 genes), gene frequency assessed, and panel performance restricted to genes in >65% of panels performed (6 genes, α = 0). Detailed descriptions of the algorithm and statistical methodology are provided in the Supplementary Data.

Analyses between two groups were performed with one-tailed Student t test or Mann-Whitney t test, and for more than two groups, the one-way ANOVA with Šidák or Bonferroni post hoc test was performed. P < 0.05 was considered significant. Error bars in figures represent the SD.

To generate Tregs and Tconvs in multiple states of activation, naive (CD4⁺CD25^hiCD45RA⁺) and memory (CD4⁺CD25^hiCD45RA⁻) Tregs and naive (CD4⁺CD25^-CD45RA⁺) and memory (CD4⁺CD25⁻CD45RA⁻) Tconvs were sorted from blood of seven healthy adults (Fig. 1A). RNA was isolated ex vivo or after stimulation for 40 h, promoting activation-induced FOXP3 in Tconvs (Fig. 1B). Tregs, but not Tconvs, were suppressive in vitro (Fig. 1C).

The gene expression profile of the eight cell subsets was analyzed by Marker Selection module and Gene Set Enrichment Analysis. Hierarchical unsupervised clustering revealed that the composite expression of 31 genes distinguished Tregs from Tconvs regardless of activation state (Fig. 1D). This differential gene signature was validated using NanoString (Fig. 1E) and separated Tregs from Tconvs by using three independent Treg sorting strategies: naive CD25^hiCD45RA⁺, memory CD25^hiCD45RA⁻, or total CD25^hiCD127^lo (Fig. 1F).

Because many of the 31 signature genes are associated with Treg function (e.g., FOXP3, CTLA-4), we hypothesized that a change in Treg signature may be detected in diseases associated with Treg dysfunction, such as T1D. Therefore, age-, sex-, and BMI-matched (Fig. 2A) (data not shown) subjects with new-onset (n = 29) or established (n = 27) T1D and pediatric control subjects (n = 24) were recruited. Subjects with T1D had stable glycemic control as indicated by comparable fed-state levels of glycated hemoglobin and daily insulin dose (Fig. 2B).

CD4⁺CD25^hiCD127^lo Tregs and CD4⁺CD25^loCD127^hi Tconvs were sorted (Fig. 2C), with comparable FOXP3 expression in all three cohorts (Fig. 2C) and demethylated TSDR in Tregs (Fig. 2D). The Treg gene signature expression was measured by NanoString. PCA of the raw log₂-transformed data showed that the composite 31-gene signature effectively discriminated between Tregs and Tconvs from pediatric samples regardless of disease state (Fig. 2E and Supplementary Fig. 1). IL1R1 was below background and excluded from further analysis, leaving a 30-gene panel.

The composite 30-gene signature was then compared among subject groups using PCA, and the algorithm performance was tested (α = 0) (Fig. 3 and Supplementary Fig. 2A–C). Control versus new-onset T1D Tregs had biomarker potential, with an AUC of 0.73 ± 0.15 compared with an AUC of 0.59 ± 0.16 for control subjects versus subjects with established T1D and 0.50 ± 0.16 for subjects with new-onset versus established T1D. No significant correlations were found in Tconvs gene expression or the ratio of Tregs/Tconvs (data not shown).

To optimize classification power of a T1D-specific algorithm, the biomarker discovery analysis was repeated with α = 0.65 (median 10 genes/panel, error rate 0.31 ± 0.12). By using this approach, we identified six genes present in >65% of panels (Fig. 4A). Considering each of these six genes individually, three were significantly downregulated (TNFRSF1B [CD120b], TMEM23 [SGMS1], and FOXP3), ANK3 and ZNF532 were not significantly changed, and one gene was slightly upregulated (LRRC32 [GARP]) in new-onset T1D compared with control Tregs (Fig. 4B). Despite the difference in FOXP3 mRNA, FOXP3 protein levels were not different (Fig. 2C), and indeed, due to multiple layers of posttranscriptional regulation, differences in mRNA expression are not necessarily linked to simultaneous differences in protein expression (20). PCA depicting expression of these six genes showed distinct clustering of new-onset T1D versus control Tregs (Fig. 4C). Tregs from established T1D had intermediate gene expression, with only FOXP3 and ANK3 expression significantly different from that of control subjects.

To explore the utility of these six genes as a biomarker algorithm in T1D, the performance test was repeated with α = 0, restricting panels to these six genes. The six-gene panel improved biomarker performance, resulting in an AUC of 0.84 ± 0.12 and sensitivity of 0.86 ± 0.14 (Fig. 4D and Supplementary Fig. 2A, D, and E), correctly classifying 26 of 29 subjects with new-onset T1D and 20 of 24 control subjects across 1,000 panels (Fig. 4E). In contrast, the gene expression from established T1D Tregs was variable but significantly different from control Tregs, possibly reflecting the different stages of autoimmune activity in these subjects.

Changes in Treg function are believed to contribute to T1D, but tools to accurately monitor Tregs throughout the disease course are lacking. We describe a novel gene signature of human Tregs and a refined six-gene algorithm as a sensitive and specific biomarker in new-onset T1D. The data show that monitoring Treg gene expression could potentially be used as a clinically applicable biomarker to track the fitness of Tregs, overcoming many of the current limitations of protein-based markers and functional assays.

Previous attempts to define specific markers of human Tregs focused on comparing CD25⁺ and CD25⁻ cells (11,12,21,22), thereby introducing bias toward identifying T-cell activation/memory markers rather than Tregs. We developed a new 31-gene signature of human Tregs that can, when used in composite, effectively discriminate between Tregs and Tconvs regardless of their activation state. Although some of the Treg-specific genes in this signature overlap with previously identified ones (21,22), some do not, illustrating the importance of integrating data from activated Tconvs. This signature has broad applicability to T1D and many other diseases because it represents the first tool to define Tregs not limited by similar expression patterns in activated Tconvs.

Although Tregs are believed to be dysfunctional in T1D, defining whether this dysfunction relates to changes in number, intrinsic function, and/or the susceptibility of other immune cells to Treg suppression has been difficult. The proportion of CD25^hiCD127^loFOXP3⁺ Tregs in blood of subjects with T1D is normal (3); we confirmed this in a pediatric cohort (including similar levels of FOXP3 protein expression) and further showed that ex vivo Tregs from subjects with T1D have a demethylated TSDR. Whether peripheral T1D Tregs are able to suppress effector cell proliferation is unclear (4,5,23). Subtle changes in function may only be revealed under specific conditions (4), in subfractions of Tregs (6,7), and/or in subjects with specific genetic risk factors (24). For example, single nucleotide polymorphisms affecting IL-2R signaling (CD25 rs12722495, rs3118470, rs2104286) or the function of PTPN2 (rs478582, rs45450798, rs1893217) are associated with changes in Treg function (24–29), so future studies to correlate changes in the Treg signature with genetic risk factors are warranted.

Traditional assays to test Treg function are not well suited for routine clinical implementation to monitor Tregs over time because they are labor intensive, variable, and may not accurately represent in vivo function. In the current study, we show that measuring the expression of a few genes in Tregs sorted from blood represents a simple and sensitive method to detect Treg changes in children with new-onset T1D. Tregs from subjects with new-onset T1D had significantly lower expression of FOXP3, TNFRSF1B (CD120b), and TMEM23 (SGMS1), a small decrease in ZNF532 and ANK3, and slightly increased LRRC32 (GARP). These changes in gene expression are consistent with the notion that functional changes in T1D Tregs are subtle and related to FOXP3 expression, signaling, and stability (3). Preliminary experiments revealed that the six-gene algorithm may also be applicable to unfractionated PBMCs. We found an AUC of 0.64 for control subjects (n = 8) versus subjects with new-onset T1D (n = 6). Future studies with a larger sample size will be needed to fully test the possibility of using the Treg gene signature in PBMCs, which would be ideal for clinical application.

Both TNFRSF1B (CD120b) and LRRC32 (GARP) have been previously shown to be associated with Treg function. There is some evidence that CD120b may mediate some of the pro-Treg functions of tumor necrosis factor-α. For example, tumor necrosis factor-α stimulates Treg proliferation through CD120b (30,31), CD120b signaling stabilizes Tregs in inflammatory environments (31), and expression of CD120b on Tregs is necessary for suppression of colitis (32). GARP is a receptor for latent transforming growth factor-β and can transfer the mature transforming growth factor-β to its receptor on the same cell or another target cell (33). GARP is upregulated upon T-cell receptor activation (34) and is part of a positive feedback loop with FOXP3 (35). Thus, three of the six signature genes have been directly linked to Treg function and stability. The other three might also be involved in Treg signaling because they have been associated with processes important for signal transduction and immunological synapse formation (TMEM23 [SGMS1]), cellular arrangement (ANK3 is integral to the cytoskeleton), and transcriptional control (ZNF532) (36).

The Treg gene signature in subjects with established T1D was variable, with ∼40% of subjects possessing a signature similar to control subjects and the other ∼60% being similar to cells from subjects with new-onset T1D. Patients with established T1D with a disease-classified Treg signature may be suited for Treg-targeted therapy, whereas patients with a control-classified Treg signature might have quiescent autoimmunity. More research in prospective longitudinal samples (e.g., from clinical trials studying the effects of new immunomodulatory treatments) is required to test the possibility that defects in self-tolerance may be more readily detected early in disease when fulminant autoimmunity exists. Additionally, testing of prediabetic samples would be of use to determine whether the Treg signature algorithm predicts progression to diabetes. Alterations in the Treg gene signature also may be detectable more broadly in autoimmunity and/or other immune-mediated diseases, with a correlation to the presence of genetic susceptibility to autoimmunity. Whether similar or distinct changes in Treg gene expression occur in other autoimmune diseases and, if so, whether they are linked to genetic risk and/or change with disease progression/treatment are of significant interest.

In conclusion, we developed a tool to discriminate between Tregs and Tconvs regardless of their activation state, which will have broad utility in monitoring and identifying human Tregs. Because changes in expression of a few genes in the Treg signature accurately classify subjects with new-onset T1D, this assay could be used as a biomarker to track Treg changes over the disease course and to identify individuals who would be best suited for Treg-targeted therapies.

Acknowledgments. The authors thank Susan Farrell (British Columbia Children’s Hospital) for subject recruitment, Rusung Tan and Huilian Qin (Child & Family Research Institute [CFRI] and The University of British Columbia) for contributions to the diabetes biobank, and Paul C. Orban and C. Bruce Verchere (CFRI and The University of British Columbia) for ongoing support and discussions.

Funding. This study was supported by grants from the Canadian Institutes of Health Research (CIHR) (MOP-115199), the Canucks for Kids Fund, Genome British Columbia, Genome Québec, and STEMCELL Technologies Inc. A.M.P. holds fellowships from the CIHR Strategic Training Initiative in Health Research Program in Transplantation and the 4 What Matters Foundation and a JDRF postdoctoral fellowship. A.Y.W. was funded by a CFRI Research Methodology Training Grant. A.S. is supported by a CIHR Banting and Best Doctoral Award. M.K.L. receives a Scientist Salary Award from CFRI.

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

Author Contributions. A.M.P. and A.Y.W. contributed to the study design, execution of the experiments, and writing of the manuscript. A.S. performed the bioinformatic and statistical analysis. J.G., Y.K., and D.N. contributed to the sample processing and analysis. C.A.P. contributed to the study design. W.N.H. supervised the gene expression analysis and contributed to the study design. S.J.T. supervised the bioinformatic and statistical analysis and contributed to the study design and writing of the manuscript. C.P. oversaw subject recruitment and sample biobanking and contributed to the writing of the manuscript. M.K.L. contributed to the overall design and execution of the research and writing of the manuscript. M.K.L. 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.