Human immune phenotyping reveals accelerated aging in type 1 diabetes

Abstract

The proportions and phenotypes of immune cell subsets in peripheral blood undergo continual and dramatic remodeling throughout the human life span, which complicates efforts to identify disease-associated immune signatures in type 1 diabetes (T1D). We conducted cross-sectional flow cytometric immune profiling on peripheral blood from 826 individuals (stage 3 T1D, their first-degree relatives, those with ≥2 islet autoantibodies, and autoantibody-negative unaffected controls). We constructed an immune age predictive model in unaffected participants and observed accelerated immune aging in T1D. We used generalized additive models for location, shape, and scale to obtain age-corrected data for flow cytometry and complete blood count readouts, which can be visualized in our interactive portal (ImmScape); 46 parameters were significantly associated with age only, 25 with T1D only, and 23 with both age and T1D. Phenotypes associated with accelerated immunological aging in T1D included increased CXCR3+ and programmed cell death 1-positive (PD-1+) frequencies in naive and memory T cell subsets, despite reduced PD-1 expression levels on memory T cells. Phenotypes associated with T1D after age correction were predictive of T1D status. Our findings demonstrate advanced immune aging in T1D and highlight disease-associated phenotypes for biomarker monitoring and therapeutic interventions.

Keywords: Autoimmune diseases; Autoimmunity; Diabetes; Immunology.

Figure 1. Immune population dynamics and QC of outcome measures.

(A) Schematic representation of hierarchical gating strategy used to identify 172 immune cell subsets evaluated from human peripheral blood. An additional 20 parameters were derived from CBC. (B) Low technical variation observed from peripheral blood samples (n = 12) stained in duplicate for assessment by flow cytometry. (C) –log10(p) (D) and correlation strength (absolute value of ρ) from Spearman’s correlation between phenotypes (each phenotype is a data point) and age showing strongest associations in the adaptive compartment as compared with innate or CBC. Kruskal-Wallis test with Dunn’s multiple-comparison test results denoted above bars. (E–L) Phenotype proportions estimated using a smoothing spline model as a function of age in AAb- individuals. (See Supplemental Figures 1–6 and Supplemental Tables 1 and 7.) QC, quality control; CBC, complete blood count; Tcm, T central memory; Tem, T effector memory; Temra, T effector memory CD45RA⁺; Tfh, T follicular helper; Tconv, T conventional; MNC, mononuclear cells; PBPC, plasmablasts/plasma cells; DN, double negative; DP, double positive.

Figure 2. Immunophenotype trajectories in T1D.

(A) Heatmap of smoothed phenotype trajectories as a function of age in AAb- individuals with analysis restricted between the ages of 5 and 75 years to avoid predicting from sparse data. The age distribution of the cohort within this age range is shown (top histogram). Immune cell phenotypes were clustered into 4 distinct groups (axis colors, right) using hierarchical clustering (dendrogram, left). (B) Line plots of each smoothed phenotype as a function of age demonstrate distinct dynamic behavior within the 4 clusters. (C) Heatmap of smoothed phenotype trajectories as a function of age in T1D individuals with the rows arranged as in A. (D) Line plots of each smoothed phenotype as in B with the T1D smoothed phenotypes overlaid in red. (See Supplemental Figures 7 and 8.) Shifts in cluster trajectories for T1D versus AAb- were compared using a 2-tailed t test (cluster 1, P < 0.001; cluster 3, P = 0.034).

Figure 3. Immunophenotype age modeling reveals accelerated aging in T1D.

(A) Averaged coefficients from the random lasso model for all phenotypes above an empirically estimated threshold, showing those increasing with age (yellow) and decreasing with age (gray). (B) The random lasso model was used to estimate immunological predicted age in CTR (gray), T1D (red), and REL (blue). The correspondence of predicted age with chronological age is shown using a piece-wise regression model with a break at chronological age 30. (C) Residual immunological age is calculated from a linear regression of predicted age and chronological age (<30 years, n = 193). Partial regression plots between residual age and (D) BMI percentile, (E) T1D duration, and (F) rested blood glucose are shown for the multivariable regression model, along with the standardized coefficient and P value (<30 years of age, n = 90). (See Supplemental Table 3, Figure 5, and Supplemental Figures 10 and 11.)

Figure 4. Age-corrected phenotypes reveal T1D-specific differences.

As an example of the utility of our model, (A) in uncorrected data, there is no significant difference detected between T1D (n = 232) and CTR (n = 240). (B) Using the GAMLSS-corrected data, there is a significant difference between T1D and CTR. (C) All age-corrected data are available for download and analysis via the ImmScape R/Shiny application. (D) Age-corrected quantile values for T1D versus CTR were compared using nonparametric Kruskal-Wallis test and post hoc Dunn’s test with Benjamini-Hochberg multiplicity adjustment. Phenotypes increased (red) and decreased (blue) in T1D (regardless of age) are shown. (See Supplemental Table 4, Figure 5, and Supplemental Figures 9–11.) MCV, mean corpuscular volume; MCH, mean corpuscular hemoglobin; RDW, red cell distribution width; PD1, programmed cell death 1.

Figure 5. Age- and T1D-associated phenotypes.

Rectangular Venn diagram summarizes phenotypes with unique association to age (left, blue shading) or T1D (right, yellow shading), with common phenotypes displayed in the overlapping area (center, green shading). The total number of phenotypes that are “unique” or “common” to age or T1D are indicated in parentheses. A color bar illustrating the magnitude and direction of effect for age or T1D is to the left of each phenotype (bar length represents the effect size; bar color indicates the phenotype is upregulated [yellow or red] or downregulated [gray or blue] in age or T1D, respectively). (See Figure 3A and Figure 4D.)

Figure 6. Increase in CXCR3-expressing T cell subsets and increased frequency, albeit lower intensity, of PD-1 expression in T1D.

Age-corrected quantile values for (A–E) CXCR3^lo, CXCR3^–, or CXCR3⁺ frequency; (F–H) PD-1⁺ frequency; and (I–L) PD-1 MFI on T cell phenotypes (CTR n = 240, REL n = 293, RSK n = 23, T1D n = 232) for A–D; (CTR n = 247, REL n = 299, RSK n =24, T1D n = 235) for E; (CTR n = 247, REL n = 298, RSK n = 23, T1D n = 237) for F–L. Significant P values shown on graph (Kruskal-Wallis test with post hoc Dunn’s test and Benjamini-Hochberg multiplicity adjustment). (See Supplemental Figure 14.)

Figure 7. Increased HLA expression on monocytes in HLA-DR4 individuals.

(A) Volcano plot showing QTL analysis results of all flow cytometry phenotypes versus T1D risk loci. Associations shown according to direction and effect size (β) of each SNP on T1D risk. Blue designates higher and black designates lower data density. Associations between HLA-DR MFI on monocytes and tag SNPs for HLA-DR4 (rs7454108), -DR3 (rs2187668), and -DR15 (rs3129889) T1D risk or protective class II HLA alleles highlighted in red. (B) The GAMLSS model fit on all AAb- (CTR and REL combined, n = 562) to correct for age. Quantiles of HLA-DR MFI on monocytes in (C) whole cohort (n = 806), (D) CTR (n = 248), (E) REL (n = 302), (F) RSK (n = 24), (G) and T1D participants (n = 232) according to number of copies of HLA-DR4. Significant P values shown on graph (Kruskal-Wallis test with post hoc Dunn’s test and Benjamini-Hochberg multiplicity adjustment). (See Supplemental Table 6.) X, any HLA-DR allele other than DR4.