Treg gene signatures predict and measure type 1 diabetes trajectory

Abstract

Introduction

Insufficient FOXP3⁺ Treg control of T cell–mediated destruction of β cells likely contributes to type 1 diabetes (T1D) (1). Monogenic mutations in FOXP3 resulting in a lack of functional Tregs lead to T1D (2), and several T1D risk alleles occur in Treg-associated genes (e.g., SNPs in CD25, PTPN2, and PTPN22) (1). Accordingly, multiple studies are testing whether immunotherapy-based methods to boost Treg function can prevent or delay T1D progression (1, 3, 4). For example, inflammatory cytokine blockade (e.g., ustekinumab, αIL-12/23p40; ref. 5), costimulation targeting (e.g., alefacept, LFA3-4-Ig; ref. 6), or Treg expansion (e.g., low-dose IL-2; cellular therapy; refs. 7, 8) aims to restore the immunoregulatory balance. However, the availability of clinically applicable tests to measure changes in immunity during such trials is limited (3); the search for biomarkers to stratify patients likely to respond to a given therapy and track changes in immune regulation is ongoing (4, 9).

We previously developed a composite biomarker assay that measures expression of 37 genes and discriminates between Tregs and conventional T cells (Tconvs) regardless of activation state and showed that Tregs from pediatric new-onset T1D patients have a significantly altered gene signature (10). Here, we tested whether Treg gene signatures (TGSs) could also identify Treg alterations in adults with T1D and investigated the predictive power of TGS-based biomarkers to predict the rate of C-peptide decline.

Results

Adults with T1D have an altered TGS in sorted Tregs and PBMCs.

We previously found that pediatric new-onset T1D Tregs had a significantly different TGS compared with healthy children (10) and tested whether this was also true in adults. Tregs (CD4⁺CD25^hiCD127^lo) were isolated from young adults with new-onset T1D (18–35 years of age, <100 days after diagnosis, n = 20, Figure 1A), a cross-sectional T1D cohort (ref. 9, n = 40, Figure 1B), or age- and sex-matched healthy controls (HCs, n = 10 or n = 37). Expression of the TGS was measured using a nCounter FLEX, and biomarker discovery analysis was applied. We found that the TGS distinguished between T1D and healthy Tregs with high area under the receiving operating characteristic curves (AUCs) in both cohorts (Figure 1, A and B, AUC = 0.830 and = 0.953, respectively). Cutoffs for sensitivity (true positive) and specificity (true negative) determination were set between 0.25 and 0.75, revealing high sensitivity and specificity for both cohorts.

Open in a separate window

Figure 1

Altered intrinsic Treg gene signature in adults with T1D.

CD25^hiCD127^lo Tregs were sorted from (A) adults with new-onset T1D (n = 20) or (B) a cross-sectional T1D cohort (n = 41) and age- and sex-matched controls (n = 10 for A or n = 37 for B). The Treg gene signature in the resulting lysates was measured on the nCounter FLEX platform. Biomarker scores and details (model of biomarker analysis selected [elastic net {en}], number [n] of genes, AUC, sensitivity, and specificity) for the best algorithm differentiating between T1D and healthy controls are indicated next to biomarker score dot plots. Horizontal lines represent means, with SD represented by error bars; dashed horizontal lines represent cutoffs for sensitivity and specificity calculations.

Genetic risk alleles, including SNPs in CD25, PTPN2, and PTPN22, are associated with CD25 expression, response to IL-2, FOXP3 stability, and Treg function in T1D (1). We thus hypothesized that genotype information might refine biomarker accuracy. Biomarker discovery using the cross-sectional T1D cohort and including genotype(s) revealed that while overall AUCs were similar to those obtained without genotype, there was slightly improved sensitivity in several cases (Figure 2). For example, adding CD25 rs2104286 number of disease variants improved the sensitivity from 0.875 (Figure 1B) to 0.950 (Figure 2C), but overall inclusion of SNPs only modestly improved accurate discrimination between healthy and T1D Tregs.

Open in a separate window

Figure 2

Addition of Treg-associated T1D SNPs enhances sensitivity of Treg-based algorithms.

New algorithms were created using combinations of the gene signature data from sorted Tregs (cross-sectional T1D and healthy cohorts), and a single indicated SNP genotype was included as a categorical value (i.e., each genotype given a specific weight) or as the number of disease variants/alleles (NDV) (i.e., 0, 1, or 2) encoding the indicated SNP. Shown are the best algorithms for each SNP (A) genotype included as categorical value and (B–H) included as NDVs. (I) Biomarker score and algorithm for gene signature data combined with the average number of (Avg) disease variants for all SNPs assessed in A–H. Horizontal lines represent means, with SD represented by error bars; dashed horizontal lines represent cutoffs for sensitivity and specificity calculations.

Clinical application of biomarker tests should ideally entail minimal sample processing to simplify implementation, reduce processing errors, and increase reproducibility. We therefore next examined whether differential TGS expression could be detected in unfractionated peripheral blood mononuclear cells (PBMCs). Gene expression in PBMC lysates from adult new-onset (n = 19, Figure 3A) and cross-sectional (n = 35, Figure 3B) T1D cohorts and age- and sex-matched controls (HC, n = 10 or n = 38) were measured, and biomarker discovery analysis was applied. Surprisingly, gene expression in PBMC lysates also revealed highly sensitive and specific algorithms that discriminated between HCs and new-onset or cross-sectional T1D cohorts (AUC = 0.895 or = 0.977, respectively). Adding genotype information to biomarker discovery generated near perfect classification of T1D versus HCs with most AUCs, sensitivity, and specificity >0.900 (Figure 4). Thus, the inclusion of Treg-associated SNPs is valuable when assessing changes in the overall immunoregulatory balance by measuring the TGS in PBMC lysates.

Open in a separate window

Figure 3

Detection of an altered Treg gene signature in PBMCs from adults with T1D.

The Treg gene signature was measured in PBMCs from (A) adults with new-onset T1D (n = 19) or (B) a cross-sectional T1D cohort (n = 36) and age- and sex-matched controls (n = 10 for A or n = 37 for B). Biomarker scores and details (as described in Figure 1; differential expression was assessed using a 2-tailed unpaired moderated t test from the LiMMa R package [ref. 22], with a P value threshold of 0.1) for the best algorithm differentiating between T1D and healthy controls are indicated next to biomarker score dot plots. Horizontal lines represent means, with SD represented by error bars; dashed horizontal lines represent cutoffs for sensitivity and specificity calculations.

Open in a separate window

Figure 4

Addition of Treg-associated T1D SNPs enhances sensitivity of PBMC-based algorithms.

New algorithms were created using combinations of gene signature data from PBMCs (isolated from cross-sectional T1D and healthy cohorts) and a single indicated SNP genotype included as categorical value (i.e., each genotype given a specific weight) or as the number of disease variants/alleles (NDV) (i.e., 0, 1, or 2). (A and B) Genotype was included as categorical value, and (C–H) genotype was included as NDVs. (I) Biomarker scores and algorithm details combining gene signature data and the average number of (Avg) disease variants across all SNPs assessed. Horizontal lines represent means, with SD represented by error bars; dashed horizontal lines represent cutoffs for sensitivity and specificity calculations. glm, general linearized model.

Differences in gene usage between Treg- and PBMC-derived algorithms.

The varying number and identity of genes in each best-performing algorithm led us to compare gene usage across algorithms derived from Tregs versus PBMCs in the cross-sectional T1D cohort (Figure 5). We found that there were 10 core genes in all of the best-performing algorithms derived from sorted Tregs (Figure 5A). These included genes associated with DNA accessibility (METTL7A, a methyltransferase), transcription (ZNF532, RBMS3), translation (RPL23A, a ribosomal protein), signal transduction (VAV3, TCR-mediated signaling; ref. 11), and cytokine receptors (CSF2RB, IL1R1, IL7R).

Open in a separate window

Figure 5

Differential gene usage in Treg- versus PBMC-derived algorithms.

Summary of the genes present in the best-performing algorithms from the cross-sectional T1D and healthy cohorts. (A) Treg-based algorithms described in Figure 1B and Figure 2. (B) PBMC-based algorithms described in Figure 3B and Figure 4. Each of the 37 mRNAs measured is listed; gray and white squares indicate genes that were or were not present in the best-performing algorithm, respectively.

Examination of the PBMC-based algorithm also revealed 10 core genes (Figure 5B), with notable similarities and differences to the Treg-derived algorithms. For example, FOXP3 (significantly lower expression in T1D PBMCs versus controls, P = 0.0034, data not shown) and TRIB1 (significantly higher expression in T1D PBMCs versus controls, P < 0.0001, data not shown) were included in all PBMC-based algorithms but not in any derived from sorted Tregs. In contrast, IL1RN and RPL23A were always present in both Treg- and PBMC-based algorithms. Hence, PBMC-based algorithms may detect changes in Treg to Tconv ratios that are eliminated when analyzing TGS of sorted Tregs.

TGS algorithms in type 1 versus type 2 diabetes.

TGS changes could be related to poor glucose homeostasis rather than a change in immune regulation. We thus measured the TGS in sorted Tregs or PBMCs from type 2 diabetes (T2D) patients and reran the age- and sex-matched T1D and HC lysates in parallel, in this case using the nCounter SPRINT system. Principle components analysis was performed separately on Treg and PBMC data to highlight differences between patients (i.e., HC vs. T1D vs. T2D) rather than reiterate differences present between sample types (i.e., Tregs vs. PBMCs). Analysis of the unweighted TGS revealed distinct clustering of T2D, T1D, and HCs for both Treg and PBMC samples (Figure 6).

Open in a separate window

Figure 6

Altered Treg gene signature in T2D compared with T1D and healthy controls.

The Treg gene signature was measured in sorted CD25^hiCD127^lo Tregs and PBMCs from the indicated age- and sex-matched cohorts using the nCounter SPRINT system. For Tregs: T2D, n = 29; T1D, n = 7; and healthy controls (HC), n = 9. For PBMCs: T2D, n = 33; T1D, n = 10; and HCs, n = 9. Shown is a principal component analysis representing expression of all 37 genes by sample group in (A) Tregs or (B) PBMCs. 95% confidence intervals are overlaid as ellipses. T1D, type 1 diabetes; T2D, type 2 diabetes.

For these experiments, since we reran T1D and HC cohort samples from Figure 1 and Figure 3 with the T2D samples on the nCounter SPRINT system, we also had the opportunity to test T1D versus HC algorithms from Figure 1B “off-the-shelf,” finding excellent replication of the AUC (AUC = 1.000) and, without changing cutoffs, perfect sensitivity (sensitivity = 1.000) but low specificity (specificity = 0.222, Figure 7A). New biomarker analysis comparing T1D to T2D, or T2D to HC, using Treg lysates also revealed algorithms that were highly sensitive and specific (AUCs > 0.9, Figure 7, B and C). Similar results were obtained when PBMC lysates were analyzed. “Off-the-shelf” application of the algorithm from Figure 3B revealed a high AUC (AUC = 0.922) and sensitivity (sensitivity = 0.900), but low specificity (specificity = 0.444, Figure 7D). As with sorted Tregs, new biomarker discovery in PBMCs also revealed distinct gene expression in T1D versus T2D and between T2D and HCs (Figure 7, E and F). Each cohort comparison was defined by distinct sets of genes (Figure 7, G and H), indicating there are intrinsic and extrinsic differences in Treg biology in people with T1D or T2D and HCs.

Open in a separate window

Figure 7

Treg gene signature algorithms differentiate type 1 or type 2 diabetes and healthy controls.

PBMC and Treg lysates from the indicated type 1 diabetes (T1D) or type 2 diabetes (T2D) cohorts were run on a nCounter SPRINT system together with age- and sex-matched healthy control (HC) samples. (A) Biomarker scores and performance when the algorithm from Figure 1B was applied to T1D and HC Treg data. (B and C) Biomarker scores and details (as described in Figure 1) of the best biomarker algorithm differentiating between Tregs from (B) T1D or (C) HCs and T2D samples. (D) Biomarker scores and performance when the algorithm from Figure 3B was applied to T1D and HC PBMC data. (E and F) Biomarker scores and details (as described in Figures 1 and ​and3)3) of the best biomarker algorithm differentiating between (E) T1D or (F) control and T2D samples in PBMCs. Horizontal lines represent means, with SD represented by error bars; dashed horizontal lines represent cutoffs for sensitivity and specificity calculations. (G and H) Summary of gene usage in each (G) Treg- or (H) PBMC-based algorithm described in A–F. Each of the 37 mRNAs measured is listed; gray and white squares indicate genes that were or were not present in the best-performing algorithm, respectively.

Application of the TGS as clinical trial monitoring tool.

A limitation of T1D disease management and clinical trials is the lack of noninvasive biomarkers to predict C-peptide trajectory and measure intervention benefit (9). Hence, we asked whether the TGS might predict disease trajectory. Together with the JDRF Biomarker Working Group and the Immune Tolerance Network, we obtained longitudinal samples from the placebo arms of the T1DAL (Type 1 Diabetes and Alefacept) and START (Study of Thymoglobulin to Arrest Type 1 Diabetes) trials (9). Subjects were divided into tertiles with slow, moderate, or rapid C-peptide decline, as determined by mixed meal tolerance tests (MMTTs), according to absolute change in C-peptide levels from baseline (<100 days after diagnosis) to 24 months later (Figure 8A).

Open in a separate window

Figure 8

The Treg gene signature as a predictive biomarker of C-peptide decline.

(A) C-peptide was quantified (2h AUC MMTT) in new-onset T1D patients in the placebo arms of the T1DAL and START clinical trials (see CONSORT flow diagrams in Supplemental Figures 1 and 2) at baseline (M0), 6 months (M6), and 24 months (M24). The absolute change in C-peptide from M0 to M24 was calculated, and subjects were divided into those with slow (n = 9), moderate (n = 8), or rapid (n = 7) decline based on terciles. (B and C). Tregs (CD4⁺CD25^hiCD127^lo) were sorted from cryopreserved PBMCs isolated from these subjects at baseline (M0), and the TGS was measured. (B) Biomarker score and details (as described in Figures 1 and ​and3)3) of the best algorithm predicting future rapid, slow, or moderate C-peptide decline. Horizontal lines represent means, with SD represented by error bars; dashed horizontal lines represent cutoffs for sensitivity and specificity calculations. (C) Expression of the indicated genes plotted by rate of C-peptide decline group, with 1-way ANOVA with Tukey’s multiple comparisons test (P ≤ 0.05, but data not significantly different).

We measured the TGS in Treg and PBMC lysates from month 0 samples and performed biomarker discovery, seeking an algorithm that could predict future rapid versus slow C-peptide decline. Although the sample size was small (n = 7 rapid vs. 9 slow decline), the TGS in Treg lysates predicted C-peptide decline with an AUC = 0.730 (Figure 8B). Subjects with a moderate rate of C-peptide decline showed an intermediate biomarker score. Furthermore, plotting normalized gene expression of month 0 Treg lysates showed trends in line with C-peptide decline (Figure 8C). Specifically, expression of all 4 genes incorporated in the algorithm either increased (C8ORF70, PMSL11, STAM) or decreased (ICA1), with progressively worse future C-peptide decline. Expression of 2 additional genes not incorporated in the algorithm, ID2 and ZBTB38, also showed trends in line with C-peptide decline. In contrast, biomarker analysis of PBMC lysates was ineffective at discriminating between these groups (data not shown), suggesting that Treg-intrinsic alterations rather than an overall change in immunoregulatory balance may be better able to predict future disease course.

Finally, we tested if the TGS could predict C-peptide decline in response to immunotherapy. Specifically, we used samples from a phase I/II safety trial of ustekinumab (αIL-12/23 p40) in adult new-onset T1D ({"type":"clinical-trial","attrs":{"text":"NCT02117765","term_id":"NCT02117765"}}NCT02117765). Ustekinumab is commonly used in psoriasis and inflammatory bowel disease (12), and its ability to block IL-12 and IL-23 might restore the abnormal balance of T helper subsets and/or IL-17–producing Tregs that we (13) and others (14, 15) previously described in T1D. C-peptide was measured by MMTT pretreatment (month 0) and at 1 and 12 months after ustekinumab treatment. Subjects were divided into slow (n = 11) or rapid (n = 5) C-peptide decline based on reported expected rates of decline (16) (Figure 9A). We first applied the unmodified algorithm from Figure 8B “off-the-shelf” to pretreatment Treg gene expression data, finding that this algorithm prospectively identified rapid versus slow C-peptide decline with AUC = 0.709; however, the specificity was low (specificity = 0.182, Figure 9B), indicating that cutoffs for sensitivity/specificity calculations need refinement.

Open in a separate window

Figure 9

The Treg gene signature as a predictive biomarker of C-peptide decline in T1D subjects treated with ustekinumab.

Adult new-onset T1D patients were treated with ustekinumab, as outlined in the CONSORT flow diagram in Supplemental Figure 3. (A) C-peptide was quantified (2h AUC MMTT) at baseline (M0), 1 month (M1), and 12 months (M12), and absolute change in C-peptide from M0 to M12 was calculated. Subjects were divided into those with slow (n = 11) or rapid (n = 5) decline based on the absolute decline at M12, with slow subjects defined as those who lost less than 0.3 pmol C-peptide/year. (B–E) The TGS was measured in sorted CD25^hiCD127^lo Tregs or PBMCs from M0 and M9 samples. (B) Test details when the algorithm from Figure 8B was applied to M0 Treg data. (C) Treg-based algorithm and biomarker scores for slow versus rapid C-peptide decline using relative TGS data (M9/M0 prior to log₂ transformation). (D) Treg- and (E) PBMC-based algorithm and biomarker scores using M9 TGS expression data. Horizontal lines in C–E represent means, with SD represented by error bars; dashed horizontal lines represent cutoffs for sensitivity and specificity calculations. (F) Summary of gene usage in C-peptide decline algorithms described in B–E.

We then asked whether the TGS changed over the course of ustekinumab therapy by dividing month 9 Treg gene expression data by pretreatment (month 0) data. Biomarker discovery revealed excellent predictive models (AUC = 0.818, sensitivity = 0.909, and specificity = 0.800) that identified subjects with slow versus rapid C-peptide decline (Figure 9C). Evidence that the algorithm improved when both pretreatment and month 9 data were included suggests that subjects with sustained versus rapidly declining C-peptide may have differing ustekinumab-driven changes in their TGS. As with the noninterventional longitudinal samples, algorithms based on month 0 PBMC lysates were ineffective (data not shown).

Finally, the TGS at month 9 could also differentiate between slow and rapid C-peptide decline (Figure 9, D and E). Interestingly, at this time point, both Treg- and PBMC-derived signatures yielded good algorithms (AUC = 0.900 and = 0.745, respectively, Figure 9, D and E). These data suggest that the TGS may detect ustekinumab-mediated changes in Treg-intrinsic gene expression and in the balance between Tregs and Tconvs upon immunotherapy. We found some gene usage overlap with algorithms shown in Figure 1–8 (Figure 9F).

Discussion

Here, we build on our previous finding that Tregs from children with new-onset T1D have an altered TGS (10) and show that this is also true in adult new-onset and cross-sectional T1D cohorts. The altered T1D TGS did not appear to be determined by changes in glucose homeostasis, as analysis of samples from subjects with T2D revealed no overlap in gene expression profiles. We also showed the potential of TGS monitoring as a biomarker of disease trajectory, possibly enabling patient stratification for immunotherapy and monitoring of therapy outcomes. TGS measurement is simple, requires very small amounts of blood, and integrates multiple aspects of Treg biology that would be difficult to quantify in individual assays.

Our biomarker discovery approach utilized a leave-one-out cross-validation approach to identify algorithms with high AUCs, sensitivity, and specificity for each comparison of interest, with the resulting algorithms typically only including a subset of the 37 quantified transcripts. Comparisons between different algorithms revealed interesting trends in transcript utilization in samples of sorted Tregs versus PBMCs. For example, analysis of gene expression in sorted Tregs revealed recurring cell surface receptors and cell signaling molecules, suggesting that T1D Tregs have intrinsic differences in their capability to interpret and transmit molecular signals. In contrast, FOXP3 was one of the most differentially expressed genes in PBMC-based algorithms, suggesting that these samples detect changes in Treg to Tconv ratios. This possibility is supported by the finding that FOXP3 was not present in algorithms derived from sorted Tregs, samples in which differential ratios between Treg to Tconv are eliminated.

Despite small sample sizes, we found that the TGS may be a novel approach to the long-standing challenge of predicting and measuring T1D disease trajectory and intervention effects on the immunoregulatory balance. Specifically, application of a single TGS algorithm in two independent cohorts was able to identify the majority of new-onset T1D subjects whose C-peptide was likely to decline rapidly. Interestingly, preexisting changes in Tregs seemed to correlate with future disease trajectory after ustekinumab treatment. These findings may indicate that a Treg-specific defect, instead of an immunoregulatory imbalance, dictates future disease course. Our data also suggest possible differential ustekinumab-mediated intrinsic effects on Tregs and the immunoregulatory balance in subjects with rapid versus slow C-peptide decline. Overall these data indicate that further testing of TGS-based algorithms in larger cohorts is warranted.

To date, only two other studies have reported prognostic biomarkers after T1D onset. Hessner and colleagues measured changes in PBMC gene expression upon culture in T1D plasma (17), finding that patients with a higher inflammatory signature at baseline had slower C-peptide decline in response to IL-1RA in the AIDA trial ({"type":"clinical-trial","attrs":{"text":"NCT00711503","term_id":"NCT00711503"}}NCT00711503). Similarly, Linsley et al. employed whole blood RNA-sequencing to find an increased T cell signature early after rituximab treatment ({"type":"clinical-trial","attrs":{"text":"NCT00279305","term_id":"NCT00279305"}}NCT00279305) in patients with rapid C-peptide decline after 1 year (18). Compared with these large “omic” based approaches, nanoString-based TGS measurement requires very few cells (<10,000) and minimal processing, making it easier to test in validation studies and more clinically feasible. As our understanding of T1D evolves from a linear to a relapsing-remitting autoimmunity model (19), it will be of interest to continue measuring the TGS prospectively over time, including before disease diagnosis, to determine if it may enable real-time evaluation of autoimmune activity.

In conclusion, these findings suggest that measuring TGSs could be a step toward a biomarker of immune status in T1D. Future application of our findings across multiple studies, together with development of a cross-platform and cross-chemistry standardization workflow, may lead to the development of “universal” algorithms that could be applied to identify rapid versus slow progressors, monitor T1D over time, and/or select subjects likely to respond to immunotherapy.

Methods

Author contributions

AMP contributed to the overall study design, execution of the experiments, and writing of the manuscript. VC performed bioinformatic and statistical analysis. JG assisted in Treg sorting and running the nanoString Sprint system. CS oversaw subject recruitment, sample biobanking, and patient data collection for the cross-sectional T1D and T2D cohorts. CS applied for and received longitudinal T1D samples and patient data from the Immune Tolerance Network (placebo arms of the START and T1DAL trials). CS oversaw Treg sorting for the cross-sectional T1D and longitudinal Immune Tolerance Network samples. AKM, AS, SC, TE, and JPD oversaw subject recruitment, sample biobanking, and patient data collection for ustekinumab trial samples, with AKM, AS, and SC as study coordinators and JPD, RT, TE, and MKL principal or coinvestigators. SJT supervised bioinformatic and statistical analysis and contributed to study design. MKL contributed to the overall design and execution of the research and writing of the manuscript and takes responsibility for the content of the article.

Supplementary Material

Acknowledgments

Version Changes

Footnotes

References