1. Introduction
Molecular biomarkers have their role besides other markers in characterizing biological age [
1], and in particular the aging immune system [
2]. The latter process, termed immunosenescence, describes the age-related deterioration of the immune system, which does not necessarily parallel chronological age [
3]. Several simple markers, based on cell type frequencies of, e.g., natural killer (NK) and T-cells, total and memory/naïve sub-populations of CD4pos and CD8pos T-cells, or CD8pos CD28neg T-cells were proposed as immunosenescence biomarkers [
4,
5,
6,
7,
8,
9,
10], while composite scores combine several markers for a comprehensive evaluation of immunological aging [
11,
12,
13]. Among those, the immune age metric IMM-AGE [
11] is considered as one of the most advanced immunosenescence biomarkers [
14,
15,
16,
17,
18]. IMM-AGE was developed from individual longitudinal profiles of composite multi-omics data on blood cell phenotypes, functional tests with stimulated cells, and gene-expression analyses. In clinical and non-clinical settings, it has demonstrated its predictive capacity concerning cardiovascular disease and mortality [
11], the risk of sepsis in trauma patients [
19], the age-related decline in cardiorespiratory fitness [
20] and work ability [
21], and the efficacy of vaccination against SARS-CoV-2 [
22]. Notably, all these application studies had to approximate the original IMM-AGE metric by a compatible set of gene expression or blood cell type data. Recently, we had established an approximation termed IMMAX (IMMune Age indeX) based on a small set of blood cell frequencies measured by flow cytometry, which enabled the estimation of IMM-AGE with reasonable accuracy in application studies [
20,
21].
Aging biomarkers, also termed aging clocks [
1], are usually constructed in relation to or as predictors of chronological age, employing a variety of methods like simple linear regression, multivariate statistical models or machine learning algorithms [
23]. For comparison to chronological age, they can be expressed as equivalent years-of-live (EYOL), defined as the chronological age of a reference exhibiting the same aging clock level [
24]. The reference value is commonly chosen as typical for the population under consideration, e.g. as the age-specific mean or median (50
th percentile), which in a simple approach could be derived by regressing the biomarker on chronological age [
11]. Notably, this is analogous to the established concept of expressing multivariate indices of the thermal environment on an equivalent temperature scale [
25]. Rescaling to commonly known units, i.e. years for biological or immunological aging similarly to using temperature concerning thermal stress, will not only facilitate the communication of research output to the public. It also allows for the categorization of individuals’ aging types (ageotypes) by defining age-specific reference intervals [
26], and for the easy calculation of the age gap as the difference in years between the biological and chronological age [
1]. Thus, when calculated with a biomarker of immune age, the age gap will represent the rate of immunological aging compared to a population reference, where positive age gap values indicate fast or accelerated aging, and negative values mark a slow or decelerated aging of an individual’s immune system independently from chronological age [
1]. This could also ameliorate methodical issues regarding the multicollinearity induced by the correlation of immunosenescence biomarkers with chronological age in the analysis and design of aging studies [
21]. However, corresponding data and schemes for deriving aging types from immunological biomarkers are lacking, specifically concerning the comprehensive metric IMMAX.
Therefore, in analogy to pediatric growth standards [
27], we established age-adjusted reference centiles covering the adult age range for the immunosenescence biomarker IMMAX by fitting generalized additive models for location, scale and shape (GAMLSS) [
28] and quantile regression models [
29,
30] to 1,605 observations. The data were pooled from the original IMM-AGE study [
11], the Dortmund Vital Study (DVS) [
20,
31], and a study on the efficacy of vaccination against SARS-CoV-2 (VAC) [
22]. With these IMMAX centiles, which were uncorrelated to chronological age by definition, we could define different levels of slow and fast aging types, rescale IMMAX as equivalent years-of-life (EYOL) and calculate immunological age gap values. We carried out methodical comparisons suggesting that the categorization into immunological aging types was independent from the different centile estimation algorithms, but heavily affected by the choice of the immunosenescence biomarker. While EYOL represents a rescaled version of the IMMAX biomarker, the age gap could be considered as a replacement for the IMMAX centiles expressed on a years-of-live scale, as both variables correlated almost perfectly and shared the independence from chronological age. When applied to preliminary longitudinal data from the ongoing Dortmund Vital Study (ClinicalTrials.gov Identifier: NCT05155397) [
31], we observed changes from the baseline examinations after 5-years follow-up in IMMAX and EYOL, which were consistent in magnitude with the follow-up period length. The non-significant changes of the IMMAX centiles and the age gaps suggest that the participants kept their position within the cohort during follow-up. These preliminary results will require future confirmatory analyses with the completed follow-up examinations.
3. Discussion
By adopting concepts from pediatric growth-chart modelling [
27], the age-adjusted centile estimation for the pooled molecular immunosenescence biomarker data enabled the characterization of the comprehensive immune age metric IMMAX in terms of fast or slow aging types covering the adult age range from 18 to 97 years. As the IMMAX centiles are uncorrelated with age by construction, they cannot be regarded as immune age biomarker, which have to correlate with age by definition [
1,
23]. Instead, the centiles represent the rate of aging of the individual immune system compared to a population reference independent of age [
1], which helps mitigating issues concerning multicollinearity in analysis and design of prospective studies involving both immunosenescence biomarkers and chronological age [
21]. In addition, the comparison of the biomarker values to a reference curve, commonly chosen as median curve P50 [
1,
11,
23], facilitated rescaling of IMMAX as equivalent years-of-live EYOL and calculating the immunological age gap [
1]. While EYOL is a rescaled version of the biomarker IMMAX, the age gap did not correlate with age, as did the IMMAX centiles, which in turn were almost perfectly correlated with the age gap. Thus, EYOL and age gap, both employing age as commonly known scale, may help communicating study outcomes concerning the biomarker IMMAX and the IMMAX centiles, respectively, to other researchers or to the public and policy makers. Notably, this is consistent with the established concept of expressing multivariate indices of the thermal environment on an equivalent temperature scale [
25].
The preliminary longitudinal data from the DVS featured the complemental use of these metrics with the observed change in IMMAX corresponding to the 5-year follow-up period when expressed as EYOL, facilitating the interpretation that, immunologically, the sample aged as expected given the length of the follow-up period. On the other hand, the non-significant changes in the values of the IMMAX centiles and age gap indicate that, with respect to immunological aging, the individuals kept their relative positions in the cohort. Summarized, this might suggest the interpretation of the observed changes as an adaptation to the aging process during the follow-up period, rather than as dysregulation of the immune system, which might be associated with individual excessive alterations contributing to the observed inter-individual variability [
33,
34].
However, due to the limited sample size these preliminary findings will require confirmatory analyses with the completed follow-up examinations, which might also aim at explaining the considerable inter-individual variability in immunological aging in relation to the personal, behavioral, environmental and work-related data gathered in the DVS [
31]. The methodological comparisons performed in our study suggest that such future analyses could rely on the simpler approach involving linear quantile regression (LQR), while the selection of the immune age biomarker will require careful consideration, where the comprehensive IMMAX metric deems preferable to single molecular immunosenescence biomarkers.
Author Contributions
Conceptualization, P.B., C.W. and E.W.; methodology, P.B.; software, P.B.; validation, C.M., P.G. and S.G.; formal analysis, P.B.; investigation, C.M.; resources, E.W. and C.W.; data curation, P.B., C.M., and P.G.; writing—original draft preparation, P.B., C.M. and C.W.; writing—review and editing, P.G., S.G. and E.W.; visualization, P.B. and C.M.; supervision, E.W. and C.W.; project administration, P.G. and S.G.; funding acquisition, E.W. and C.W. All authors have read and agreed to the published version of the manuscript.
Figure 1.
Moving from immune age biomarkers to aging type and age gap definition. (
a) Approximation of the comprehensive IMM-AGE metric [
11] by the IMMune Age indeX IMMAX [
20], with dashed line of identity, Pearson correlation coefficient (R) and root-mean-squared prediction error (RMSE). (
b) IMMAX related to age with regression lines in the IMM-AGE study, Dortmund Vital Study (DVS) and vaccination study (VAC), respectively, and with the age distributions from the three study groups overlaid on top. (
c) IMMAX distribution conditional on age with the line connecting the median values estimated for the pooled training data (N=1,580) by a generalized additive model for location, scale and shape (GAMLSS). (
d) IMMAX centiles estimated by GAMLSS in relation to age and IMMAX, respectively. (
e) Percentile curves (PXX, XX denoting the percentile) estimated by GAMLSS for the rating of IMMAX values into six aging type categories and for re-scaling IMMAX to equivalent years of life (EYOL) with the median curve (P50) as reference, and the immunological age gap defined by the difference between EYOL and chronological age. (
f) Spearman rank correlations (ρ) of IMMAX and IMMAX centiles with EYOL and age gap, respectively.
Figure 1.
Moving from immune age biomarkers to aging type and age gap definition. (
a) Approximation of the comprehensive IMM-AGE metric [
11] by the IMMune Age indeX IMMAX [
20], with dashed line of identity, Pearson correlation coefficient (R) and root-mean-squared prediction error (RMSE). (
b) IMMAX related to age with regression lines in the IMM-AGE study, Dortmund Vital Study (DVS) and vaccination study (VAC), respectively, and with the age distributions from the three study groups overlaid on top. (
c) IMMAX distribution conditional on age with the line connecting the median values estimated for the pooled training data (N=1,580) by a generalized additive model for location, scale and shape (GAMLSS). (
d) IMMAX centiles estimated by GAMLSS in relation to age and IMMAX, respectively. (
e) Percentile curves (PXX, XX denoting the percentile) estimated by GAMLSS for the rating of IMMAX values into six aging type categories and for re-scaling IMMAX to equivalent years of life (EYOL) with the median curve (P50) as reference, and the immunological age gap defined by the difference between EYOL and chronological age. (
f) Spearman rank correlations (ρ) of IMMAX and IMMAX centiles with EYOL and age gap, respectively.
Figure 2.
Influence of centile estimation models on the rating of immune aging types. (a) Age-depending percentile curves (PXX, with XX denoting the percentile) for IMMAX with corresponding aging types estimated by GAMLSS, linear (LQR) and non-parametric quantile regression (NQR), respectively. (b) Confusion matrices with the proportion (prop) of immune aging type combinations rated by GAMLSS compared to the quantile regression models LQR and NQR, respectively, for the training (upper panel) and test data (lower panel). (c) Krippendorff’s α with error bars indicating 95%-CI for assessing the agreement between the immune aging types rated by pairwise and triple combinations of the different models for both training and test data.
Figure 2.
Influence of centile estimation models on the rating of immune aging types. (a) Age-depending percentile curves (PXX, with XX denoting the percentile) for IMMAX with corresponding aging types estimated by GAMLSS, linear (LQR) and non-parametric quantile regression (NQR), respectively. (b) Confusion matrices with the proportion (prop) of immune aging type combinations rated by GAMLSS compared to the quantile regression models LQR and NQR, respectively, for the training (upper panel) and test data (lower panel). (c) Krippendorff’s α with error bars indicating 95%-CI for assessing the agreement between the immune aging types rated by pairwise and triple combinations of the different models for both training and test data.
Figure 3.
Influence of centile estimation models and choice of biomarker on immune aging types. (a) Age-depending percentile curves (PXX, with XX denoting the percentile) estimated by GAMLSS, linear (LQR) and non-parametric quantile regression (NQR) for the five immunosenescence biomarkers predicting IMMAX: the cell frequency ratios of memory to naïve CD4 and CD8 T-cells, respectively, of NK- to T-cells, of CD4pos to CD8pos T-cells, and the frequency of CD28neg-CD8pos T-cells. (b) Bivariate Spearman rank correlations (ρ) between the immune aging types rated by 18 algorithms R1-R18 (6 biomarkers × 3 models) for the training data (upper triangular matrix) and test data (lower triangular matrix). (c) Krippendorff’s α assessing the agreement between the immune aging type rating by the different models for the six biomarkers (left panel), as well as by the different biomarkers in relation to the three models (right panel), for both training and test data.
Figure 3.
Influence of centile estimation models and choice of biomarker on immune aging types. (a) Age-depending percentile curves (PXX, with XX denoting the percentile) estimated by GAMLSS, linear (LQR) and non-parametric quantile regression (NQR) for the five immunosenescence biomarkers predicting IMMAX: the cell frequency ratios of memory to naïve CD4 and CD8 T-cells, respectively, of NK- to T-cells, of CD4pos to CD8pos T-cells, and the frequency of CD28neg-CD8pos T-cells. (b) Bivariate Spearman rank correlations (ρ) between the immune aging types rated by 18 algorithms R1-R18 (6 biomarkers × 3 models) for the training data (upper triangular matrix) and test data (lower triangular matrix). (c) Krippendorff’s α assessing the agreement between the immune aging type rating by the different models for the six biomarkers (left panel), as well as by the different biomarkers in relation to the three models (right panel), for both training and test data.
Figure 4.
Influence of reference curves, defined as 50th percentile (P50) estimated by GAMLSS and LQR, respectively, or assuming a hypothetical reference with percentiles gradually progressing with age (GAMLSS-HYP), on immunological age gap. Spearman rank correlations (ρ) for training and test data (a) between age gaps defined by three reference curves, and (b) of age gap vs. chronological age with overlaid smoothing splines.
Figure 4.
Influence of reference curves, defined as 50th percentile (P50) estimated by GAMLSS and LQR, respectively, or assuming a hypothetical reference with percentiles gradually progressing with age (GAMLSS-HYP), on immunological age gap. Spearman rank correlations (ρ) for training and test data (a) between age gaps defined by three reference curves, and (b) of age gap vs. chronological age with overlaid smoothing splines.
Figure 5.
Preliminary longitudinal data from the DVS. (a) Baseline and follow-up values from 53 participants of IMMAX (upper panel) and IMMAX centiles (lower panel) with statistical significance assessed by the Wilcoxon paired sample test. (b) Change from baseline of equivalent years of life (ΔEYOL) derived from IMMAX with different reference curves defined as 50th percentile estimated by GAMLSS and LQR, respectively, or assuming a hypothetical reference with percentiles gradually progressing with age (GAMLSS-HYP). Triangles and error bars indicate means with 95%-CI, the dashed horizontal line marks the 5 years follow-up period, corresponding to the zero-value for the change in age gap (ΔAge gap). Pairwise statistical comparison by Wilcoxon paired sample test adjusted for multiple testing. (c) Pairwise scatterplots with Spearman correlation coefficients (ρ) for the change from baseline (Δ) in IMMAX, IMMAX centile, and EYOL calculated with the GAMLSS-P50 reference; (d) Density plots comparing the female and male distributions of ΔIMMAX (left panel, with upper horizontal axis showing the linearly rescaled equivalent years of life ΔEYOLLQR-P50), and ΔIMMAX centiles (right panel). Vertical dot-dashed lines indicate medians by sex; solid lines mark zero effects. Statistical comparisons by independent-sample Wilcoxon test.
Figure 5.
Preliminary longitudinal data from the DVS. (a) Baseline and follow-up values from 53 participants of IMMAX (upper panel) and IMMAX centiles (lower panel) with statistical significance assessed by the Wilcoxon paired sample test. (b) Change from baseline of equivalent years of life (ΔEYOL) derived from IMMAX with different reference curves defined as 50th percentile estimated by GAMLSS and LQR, respectively, or assuming a hypothetical reference with percentiles gradually progressing with age (GAMLSS-HYP). Triangles and error bars indicate means with 95%-CI, the dashed horizontal line marks the 5 years follow-up period, corresponding to the zero-value for the change in age gap (ΔAge gap). Pairwise statistical comparison by Wilcoxon paired sample test adjusted for multiple testing. (c) Pairwise scatterplots with Spearman correlation coefficients (ρ) for the change from baseline (Δ) in IMMAX, IMMAX centile, and EYOL calculated with the GAMLSS-P50 reference; (d) Density plots comparing the female and male distributions of ΔIMMAX (left panel, with upper horizontal axis showing the linearly rescaled equivalent years of life ΔEYOLLQR-P50), and ΔIMMAX centiles (right panel). Vertical dot-dashed lines indicate medians by sex; solid lines mark zero effects. Statistical comparisons by independent-sample Wilcoxon test.
Figure 6.
Gating strategy for derivation of immunosenescence biomarkers by flow cytometry. EDTA blood was stained with a Fixable Viability Dye (FVD) and for the indicated markers, which were calculated from relative cell frequencies: ratios of NK- to T-cells, of CD4pos to CD8pos T-cells, of memory to naïve CD4 and CD8 T-cells, respectively, and the frequency of CD28neg-CD8pos T-cells.
Figure 6.
Gating strategy for derivation of immunosenescence biomarkers by flow cytometry. EDTA blood was stained with a Fixable Viability Dye (FVD) and for the indicated markers, which were calculated from relative cell frequencies: ratios of NK- to T-cells, of CD4pos to CD8pos T-cells, of memory to naïve CD4 and CD8 T-cells, respectively, and the frequency of CD28neg-CD8pos T-cells.
Table 1.
Materials for PBMC staining, including antigens, antibody clones and coupled fluorochromes, distributors and antibody dilution used to stain 0.2 × 106 PBMC (DVS) or 100 µl whole blood (VAC).
Table 1.
Materials for PBMC staining, including antigens, antibody clones and coupled fluorochromes, distributors and antibody dilution used to stain 0.2 × 106 PBMC (DVS) or 100 µl whole blood (VAC).
Sample |
Antigen |
Clone |
Fluorochrome |
Company |
Dilution 1/x |
DVS |
live / dead |
|
zombie Yellow |
Biolegend |
1000 |
|
CD3 |
UCHT1 |
BV510 |
BD Horizon™ |
400 |
|
CD56 |
B159 |
PE-CF594 |
BD Pharmingen™ |
100 |
|
CD4 |
RPA-T4 |
APC-H7 |
BD Pharmingen™ |
100 |
|
CD8 |
RPA-T8 |
FITC |
BD Pharmingen™ |
200 |
|
CD197 (CCR7) |
150503 |
Alexa Fluor® 647 |
BD Pharmingen™ |
50 |
|
CD45RA |
HI100 |
Alexa Fluor® 700 |
BD Pharmingen™ |
400 |
|
CD28 |
CD28.2 |
PerCP-Cy™5.5 |
BD Pharmingen™ |
100 |
VAC |
live / dead |
|
Fixable Viability Dye eFluor™ 780 |
ThermoFisher Scientific |
400 |
|
CD3 |
UCHT1 |
BV510 |
BD Horizon™ |
100 |
|
CD56 |
B159 |
PE-Cy™5 |
BD Pharmingen™ |
50 |
|
CD4 |
RPA-T4 |
BV421 |
BD Horizon™ |
100 |
|
CD8 |
RPA-T8 |
BB515 |
BD Horizon™ |
400 |
|
CD197 (CCR7) |
3D12 |
PE |
BD Pharmingen™ |
100 |
|
CD45RA |
HI100 |
Alexa Fluor® 700 |
BD Pharmingen™ |
100 |
|
CD28 |
CD28.2 |
PerCP-Cy™5.5 |
BD Pharmingen™ |
100 |