Fig. 1 | Overview of single-cell transcriptomes of peripheral immune cells across each stage of the human lifespan.

a. Workflow of the overall experimental design of this study.

b. UMAP plot illustrating the peripheral immune cell subsets identified within PBMCs. The 25 immune cell subsets are labeled with different colors.

c. Heatmaps showing the row-scaled expression of signature genes in each cell subset. Heatmap colors from blue to red represent the relative expression levels of genes. Color scale: red, high expression; blue, low expression.

d. Dot plot showing the age-related variations in proportion of the 25 immune cell subsets, across 13 age groups. Kruskal-Wallis test was used to identify cell subsets with significant difference across age groups. FDR correction was applied to account for multiple testing and the values were shown above the dot plot, an FDR < 0.05 was considered significant.

FDR: false discovery rate; mDCs, myeloid dendritic cells; PBMCs, peripheral blood mononuclear cells; pDCs, plasmacytoid dendritic cells; UMAP, uniform manifold approximation and projection.

Fig. 2 | Global transcriptome dynamics of peripheral immune cells throughout the human lifespan.

a. Bar plot showing the number of DEGs derived from scRNA-seq analyses of all cell subsets. The dashed line indicates the threshold for the number of DEGs ≥ 100.

b. Heatmap showing the functional enrichment analyses of DEGs in all cell subsets. Heatmap colors indicate the significance of the biological processes.

c. Bar plots showing the frequencies of upregulated (red) and downregulated (blue) M-DEGs derived from scRNA-seq analyses observed in all cell subsets. All M-DEGs with frequency ≥ 5 are displayed.

d. Line plots showing the expression trajectories of upregulated (left panel) and downregulated (right panel) age-related genes identified via bulk RNA-seq analyses. Data were shown in terms of gene expression z-score. Each gray line represents expression trajectory of a gene. The red/blue line represents the average expression trajectory.

e. Venn diagrams showing the number of overlapping genes between common M-DEGs identified in scRNA-seq data and age-related genes identified in bulk RNA-seq data (left: upregulated genes, right: downregulated genes). Among the 56 common M-DEG genes defined by scRNA-seq, 13 intersected with the age-related genes defined by bulk RNA-seq.

f. Heatmap showing the expression of SELL gene in seven cell subsets via scRNA-seq and CD62L protein (encoded by SELL) in corresponding cell subsets via CyTOF among eight life stages.

CyTOF, cytometry by time-of-flight; DEGs, differentially expressed genes; M-DEGs, monotonically differentially expressed genes; scRNA-seq, single-cell RNA sequencing.

Fig. 3 | Rewiring of CCIs among peripheral immune cells across the human lifespan.

a. Heatmap showing that cell-cell interaction (CCI) probability varied significantly across the lifespan. The probability for each CCI was row-scaled with a z-score. CCIs were clustered into four clusters by hierarchical clustering, and these clusters were enriched in different age groups. The significantly varied CCIs among eight different human life stages were determined by comparing the CCI probability of all cell subset pairs in each sample using the Kruskal-Wallis test. P values were adjusted with FDR, and the adjusted p value < 0.0001 was considered statistically significant.

b. Barplot showing T cells experienced the most intensive rewiring in cell–cell interactions with age.

c–f. Network graphs showing the four CCI clusters from (a). The CCIs enriched in each cluster are highlighted in red. The hub cell subset of each cluster is highlighted with larger node size and red label.

g. Network graph showing the CCI network mediated by differentially expressed ligands and receptors. The sub-network circled in red shows the immune checkpoints and the corresponding affected cell subsets.

h. Bar plot showing the top 10 biological processes enriched in the differentially expressed ligands and receptors.

i. Boxplot showing the average shortest path length and degree of immune checkpoints, and non-immune checkpoints in (g). P values were determined by a two-sided unpaired Wilcoxon test.

j. Bar plot showing, for each immune checkpoint ligand and receptor, the number of cell subsets in which they were detected as differentially expressed.

k. Heatmaps showing the scaled expression of co-stimulatory checkpoint CD27 and ICOS, and co-inhibitory checkpoint PDCD1 genes across the lifespan in the indicated cell subsets, detected by scRNA-seq.

l. Heatmaps showing the scaled expression of co-stimulatory checkpoint CD27 and ICOS, and co-inhibitory checkpoint PD1 proteins across the lifespan in the indicated cell subsets, detected by CyTOF.

CCIs, cell-cell interactions; CyTOF, cytometry by time-of-flight; scRNA-seq, single-cell RNA sequencing.

Fig. 4 | Clonal expansion of T cells across the entire lifespan.

a. Stacked bar graph showing the distribution of T-cell clone sizes in the 13 age groups.

b. Heatmap showing the row-scaled usage of significantly differentially expressed V-J genes from TCR α-chains and TCR β-chains across six high-TCE groups (2 y, 6 y, 12 y, 70 y, 80 y, and 90 y). The six age groups were clustered into two distinct groups: 1) the 2–12 y group (early-life high-TCE group) and 2) 70–90 y group (late-life high-TCE group). Red, high expression; blue, low expression.

c. Bar plot showing the distribution of T-cell clone sizes in each αβ T-cell subset.

d, e. Boxplot showing the clonal expansion levels of total αβ T-cell subsets (d) and CD8_TEM_GNLY cells (e) across 13 age groups quantified by STARTRAC-expa metric. The black line indicates LOESS regression. Spearman correlation coefficient ρ and p values are indicated.

f. Volcano plot showing the DEGs in high-TCE CD8_TEM_GNLY cells (clone size ≥ 50 cells) in the late-life high-TCE group (70 y, 80 y, and 90 y) compared with those in the early-life high-TCE group (2 y, 6 y, and 12 y). The representative genes are labeled. Red, upregulated in the late-life high-TCE group [Log2 (fold change) ≥ 0.25, adjusted p value < 0.05]; blue, upregulated in the early-life high-TCE group [Log2 (fold change) ≤ -0.25, adjusted p value < 0.05].

g. Bar plots showing the representative enriched biological processes of DEGs upregulated in high-TCE CD8_TEM_GNLY cells in early-life (left panel) and in late-life (right panel).

h. Boxplots showing the gene module scores for cytotoxicity, oxidative phosphorylation, positive regulation of mononuclear cell proliferation, and response to type II interferon, in high-TCE CD8_TEM_GNLY cells in the early-life and late-life high-TCE groups. P values were determined by Wilcoxon test.

i. Rank for TF regulons in high-TCE CD8_TEM_GNLY cells in the early-life (left panel) and late-life (right panel) high-TCE groups using the regulon specificity score. The top ten TF regulons in each group are listed and marked in blue.

DEGs, differentially expressed genes; TCE, T-cell clonal expansion; TF, transcription factor.

Fig. 5 | Aging features of naive T cells and CD8+ MAIT cells.

a. Scatter plot illustrating dynamic changes in the proportion of CD4+ naive T cells (CD4_Naive_CCR7 cells; left panel) and CD8+ naive T cells (CD8_Naive_LEF1 cells; right panel) in total αβ T cells with aging based on scRNA-seq. The black line indicates LOESS regression.

b. Scatter plot illustrating dynamic changes in the proportion of CD4+ naive T cells (CD4_Naive cells; left panel) and CD8+ naive T cells (CD8_Naive cells; right panel) in total αβ T cells with aging based on CyTOF data. The black line indicates LOESS regression.

c. Heatmap showing the row-scaled expression of M-DEGs in CD4_Naive_CCR7 cells across eight human life stages (left panel). Heatmap colors indicate relative expression levels of genes from low (blue) to high (red). Line plots show the expression trajectories of upregulated M-DEGs in CD4_Naive_CCR7 cells (right panel). Data are shown in terms of gene expression z-score. Each gray line represents the expression trajectory of a gene. The blue line represents the average expression trajectory.

d. Heatmap showing the row-scaled expression of M-DEGs in CD8_Naive_LEF1 cells across eight human life stages (left panel). Line plots showing the expression trajectories of upregulated M-DEGs in CD8_Naive_LEF1 cells (right panel). Data are shown in terms of gene expression z-score. Each gray line represents the expression trajectory of a gene. The blue line represents the average expression trajectory.

e. Dot plot showing the representative enriched biological processes of upregulated M-DEGs in CD4+ naive T cells and CD8+ naive T cells. Dot colors represent the significance of biological processes, and dot sizes indicate the gene number annotated to each biological process.

f, g. Scatter plots showing the diversity of TCR repertoires in CD4+ naive T cells (CD4_Naive_CCR7 cells; left panel) and CD8+ naive T cells (CD8_Naive_LEF1 cells; right panel) across eight human life stages. The diversity of TCR repertoires was measured by Shannon’s Entropy (f) and Inverse Simpson Index (g). The black line indicates LOESS regression.

h, i. Scatter plot illustrating the dynamic proportion changes of CD8_MAIT_SLC4A10 cells in total αβ T cells with aging based on scRNA-seq (h) and CyTOF (i) analyses across eight life stages. The black line indicates LOESS regression.

j. Scatter plot showing TCR repertoire diversity in CD8_MAIT_SLC4A10 cells across eight life stages. The diversity of TCR repertoires was measured by Shannon’s Entropy. The black line indicates LOESS regression.

k. Heatmap showing the expression of all the DEGs (left panel) and zoom-in on upregulated genes (right panel) in CD8_MAIT_SLC4A10 cells in adolescents (12–18 y) compared with those in other life periods including children (1–8 y), adults (30–60 y), and the elderly (70–90 y).

l. Bar plots showing the representative enriched biological processes of upregulated (top panel) and downregulated (bottom panel) DEGs in CD8_MAIT_SLC4A10 cells in adolescents (12–18 y) compared with those in other life periods including children (1–6 y), adults (30–60 y), and the elderly (70–90 y).

m. GSEA enrichment plot for gene set “tuberculosis” in CD8_MAIT_SLC4A10 cells in the adolescent group compared with that in other age groups. A p value < 0.05 was considered statistically significant.

n. Schematic diagram showing the design of MAIT cell’s bactericidal capability assay. Bacterial loads in lysates of infected cells after co-culture with MAIT cells from healthy individuals aged 18 years (n = 9) and 50 years (n = 9) (left panel) were measured. Bar plot showing the bacteria load after MAIT cell co-culture, shown with mean ± s.e.m. and p value determined by the two-sided unpaired t test (right panel).

CFU, colony forming unit; CyTOF, cytometry by time-of-flight; DEGs, differentially expressed genes; GSEA, gene set enrichment analysis; MAIT, mucosal-associated invariant T; M-DEGs, monotonically differentially expressed genes; scRNA-seq, single-cell RNA sequencing.

Fig. 6 | Function of the newly identified cytotoxic B-cell subset (B_BCR+GNLY+ cells).

a. Heatmap showing the functional enrichment analyses of signature genes in each B-cell subset based on the scRNA-seq analyses. The top 10 pathways in each cell subset (ranked by p value) are listed. Heatmap colors represented the significance of the biological processes.

b. Immunofluorescence of B-cell markers (IgM, IgA, and J chain) and the cytotoxic protein granulysin (GNLY) in PBMCs collected from six healthy 6-year-olds. Green represents IgM, IgA, or J chain; red represents GNLY; blue represents nuclei; and yellow indicates double-positive cells (IgM and GNLY, IgA and GNLY, or J chain and GNLY). The experiment was repeated at least three times with similar results.

c. Bar plot showing the percentage of B_BCR+GNLY+ cells in PBMCs from three infants with longitudinal follow-up, examined at one-year-old and at two-year-old via scRNA-seq analyses. R1, R2, and R3 represent different individuals.

d. Flow cytometry scatter plots showing proportions of IgM+GNLY+ cells in children (4–9 y), adults (24–30 y), and the elderly (64–72 y), without stimulation (control, top left), and day 1 after in vitro stimulation with an expansion mixture consisting of IL-21 (10 ng/ml), IL-2 (50 IU/ml), CD40L (50 ng/ml), anti-BCR (5 mg/ml) and CpG ODN 2006 (1 mg/ml) (bottom left) (n = 5 per group). Bar plots showing the percentage of IgM+GNLY+ cells (defined by flow cytometry) in the same children, adults, and elderly groups (right panel). Bar plots are shown with the mean ± s.e.m. P values were determined by one-way ANOVA followed by the Tukey post hoc test. A p value < 0.05 was considered statistically significant.

e. Representative FACS scatter plots showing the proportion of IgM+GNLY+ cells in cultured human B cells (1.5 × 106 cells/ml) on day 1 and 3 after in vitro stimulation with an expansion mixture consisting of IL-21 (10 ng/ml), IL-2 (50 IU/ml), CD40L (50 ng/ml), anti-BCR (5 mg/ml) and CpG ODN 2006 (1 mg/ml) (left panel) (n = 4 per group). Bar plots showing the percentage of IgM+GNLY+ cells in B cells from each group (right panel). Bar plot shows the mean ± s.e.m and p values were determined by the two-sided unpaired Wilcoxon test. A p value < 0.05 was considered statistically significant.

f. Bar plots showing GNLY concentration in culture supernatants of B cells on days 1 and 3 of culture (1.5 × 106 cells/ml) after in vitro stimulation with phosphate buffered saline or an expansion mixture. The concentrations of GNLY were determined by ELISA (n = 4 per group). P values were determined by the two-sided unpaired Wilcoxon test. A p value < 0.05 was considered statistically significant.

g. Schematic diagram detailing the assay on bactericidal effects against S. aureus or S. pneumoniae for isolated B cells that overexpressing GNLY by electroporation.

h. Bar plots showing the bactericidal effects of GNLY-overexpressing B cells against S. aureus (left panel) or S. pneumoniae (right panel). Error bars represent mean ± s.e.m. P values were determined by one-way ANOVA followed by the Tukey post hoc test. A p value < 0.05 was considered statistically significant.

i. Schematic diagram of gene editing strategy for transgenic mice overexpressing human GNLY (top panel). To generate B-cell-specific GNLY-overexpression (GNLY-TG) mice, R26-GNLY (NTG) mice were crossed with the mice carrying the CD19-Cre transgene. PCR analysis was performed for the confirmation of successful gene editing (bottom panel). The 195 bp band represents the GNLY gene, which was produced by the insertion of GNLY cDNA at Rosa26 locus. The 200 bp and 280 bp bands represent CD19-Cre+ and CD19-Cre-, respectively.

j. Flow cytometry scatter plots showing proportions of CD19+GNLY+ cells in PBMCs isolated from GNLY-TG and NTG mice (left and middle panels). Bar plots showing the percentage of CD19+GNLY+ cells in GNLY-TG and NTG mice (right panel). Bar plots are shown with the mean ± s.e.m. P values were determined by the two-sided unpaired Wilcoxon test. A p value < 0.05 was considered statistically significant.

k. The bactericidal effects of B cells derived from GNLY-TG and NTG mice against S. aureus or S. pneumoniae. Error bars represent mean ± s.e.m. P values were determined by the two-sided unpaired t test. A p value < 0.05 was considered statistically significant.

l. Heatmap showing the correlation of gene expression among different B-cell subsets.

m. Heatmap showing the UMI percentage of immunoglobulin heavy-chain genes (IGHA1, IGHA2, IGHD, IGHG1, IGHG2, IGHG3, IGHG4, and IGHM) in each B-cell subset. Numbers represent the actual percentage values.

n. The order of B-cell subsets along pseudotime determined by Monocle2. Each dot represents a single cell. Color gradients represent the pseudo temporal order (top panel) and different subsets (bottom panel).

o. Network plot showing the ratio of shared BCR clones between any two cell subsets. The width of the lines is proportional to the ratio.

BCR, B-cell receptor; ELISA, enzyme-linked immunosorbent assay; PBMCs, peripheral blood mononuclear cells; UMI, unique molecular identifier.

Fig. 7 | Construction of the single-cell immune age prediction model based on lifecycle-wide single-cell data.

a. Workflow of the single-cell immune age prediction model construction.

b. Dot plot showing the predictive performance of the cell type-specific immune age prediction models based on our lifecycle-wide single-cell data. The top ten cell types with the highest scores are shown. The integrated prediction model based on all 25 cell-types had the highest predictive score and the best performance.

c. Scatter plot showing the correlation of single-cell immune age (siAge) with chronological age (cAge) in the training set. Pearson correlation coefficient (PCC) and p values are indicated. The blue line indicates linear regression. The gray shadow covers the 95% confidence interval. A p value < 0.05 was considered statistically significant.

d. Scatter plot showing the correlation of siAge with cAge in the external healthy validation cohort.

e. Scatter plot showing the correlation of siAge with cAge in the external validation cohort of diseased patients with disturbed immune functions.

f. Box plot comparing the fit index distribution calculated by siAge model in the healthy validation cohort with those of the diseased validation cohorts, which include nasopharyngeal carcinoma (NPC), systemic lupus erythematosus (SLE), and Kawasaki diseases (KD) cohorts. P values were determined by the two-sided unpaired Wilcoxon test. A p value < 0.05 was considered statistically significant. The fit index was calculated as log2(siAge /cAge).

g. Bar plot showing the top 21 key genes identified by Random Forests regression. The genes are ranked in descending order of importance with respect to the accuracy of the model. The insert represents a 10-fold cross-validation error as a function of the number of input features used to regress against the chronological ages.

h. Heatmap showing the row-scaled expression of the top 21 key genes across age in corresponding cell subsets. Red, high expression; blue, low expression.

i. Bar plot showing functional enrichment of the top 21 key genes.

Extended Data Fig. 1 | Basic characteristics of the integrated dataset and signature genes of 25 peripheral immune cell subsets covering 13 age groups.

a. UMAP plots showing the 25 peripheral immune cell subsets covering 13 age groups across the entire lifespan. Cell subsets are indicated in different colors. B_BCR+GNLY+ cells are circled with red dashed lines.

b. UMAP plots showing the distribution of expression of selected signature genes for 25 cell subsets. Expression levels are color-coded, and the legends are labeled according to the log scale. Color bars from gray to blue indicate the expression levels of genes. B_BCR+GNLY+ cells are circled with red dashed lines.

c. UMAP plots showing the distribution of expression of selected signature genes for 12 T-cell subsets. Expression levels are color-coded: gray, no expression; blue, high expression.

d. Violin plots showing the distribution of expression of selected signature genes for 25 cell subsets. The columns represent the 25 cell subsets and the rows represent the selected signature genes.

e. Violin plots showing the distribution of UMI counts (left panel) and gene counts (right panel) per cell type.

UMAP, uniform manifold approximation and projection; UMI, unique molecular identifier.

Extended Data Fig. 2 | Validation of cell annotation derived from scRNA-seq data by CyTOF, Azimuth, and GSEA.

a. UMAP projection of integrated CyTOF data of peripheral immune cell subsets identified within the PBMCs, with cells colored by major immune cell types.

b. Heatmap showing the relative mean cluster gene and protein expressions in major immune cell types from scRNA-seq (left panel) and CyTOF (right panel) analyses, respectively.

c. Heatmap showing the correlation between scRNA-seq-derived major cell types and CyTOF-derived major cell types, computed from expression of scRNA-seq variable genes and corresponding CyTOF protein markers.

d. UMAP projection of T-cell subsets identified in the integrated CyTOF data, with cells colored by T-cell subsets.

e. Heatmap showing the relative mean cluster gene and protein expressions in T-cell subsets from scRNA-seq (left panel) and CyTOF (right panel) analyses, respectively.

f. Heatmap showing the correlation between the scRNA-seq-derived T-cell subsets and the CyTOF-derived T-cell subsets, computed from expression of scRNA-seq variable genes and corresponding CyTOF protein markers.

g. Heatmap showing strong consistency of our cell annotations based on scRNA-seq and cell annotations predicted by the automated annotation tool, Azimuth. Color bars from blue to red indicate the correlation.

h. GSEA plots showing immunological signature gene sets enriched by the marker genes of each T-cell subset, illustrating that our categorization of T cells based on scRNA-seq aligned with established findings.

CyTOF, cytometry by time-of-flight; GSEA, gene set enrichment analysis; scRNA-seq, single-cell RNA sequencing; UMAP, uniform manifold approximation and projection.

Extended Data Fig. 3 | Cell composition dynamics of 25 peripheral immune cell subsets across 13 age groups.

a. Stacked bar graph showing the average proportion of each cell subset within total PBMCs based on scRNA-seq analyses across 13 age groups.

b–d. Scatter plots showing the correlation between age and the proportion of each cell subset based on the scRNA-seq analyses in Fig. 1d. Spearman correlation coefficient ρ and p values are indicated. A p value < 0.05 was considered statistically significant. (b) Cell subsets significantly positively correlated with age, (c) cell subsets significantly negatively correlated with age, (d) cell subsets enriched in certain age groups.

e–g. Scatter plots showing the correlation between age and the proportion of each cell subset based on the CyTOF analyses. Spearman correlation coefficient ρ and p values are indicated. A p value < 0.05 was considered statistically significant. (e) Cell subsets significantly positively correlated with age, (f) cell subsets significantly negatively correlated with age, (g) cell subsets enriched in certain age groups.

CyTOF, cytometry by time-of-flight; PBMCs, peripheral blood mononuclear cells; scRNA-seq, single-cell RNA sequencing.

Extended Data Fig. 4 | Functional enrichment analyses of M-DEGs in peripheral immune cell subsets.

a. Bar plot showing the numbers of total M-DEGs defined by scRNA-seq in each cell subset.

b. Bar plot showing the numbers of upregulated (left panel) and downregulated (right panel) M-DEGs in each cell subset.

c. Heatmap showing the functional enrichment analyses of upregulated M-DEGs in different peripheral immune cell subsets.

d. Heatmap showing the functional enrichment analyses of downregulated M-DEGs in different peripheral immune cell subsets.

e. Functional enrichment network of the common M-DEGs in Fig. 2c. The green nodes denote biological processes; the orange nodes denote the M-DEGs of the biological processes.

f. Functional enrichment network of the common M-DEGs in Fig. 2d. The green nodes denote biological processes; the orange nodes denote the M-DEGs of the biological processes.

g. The positive correlation of LGALS1 expression with age and the negative correlation of CCR7 and TMIGD2 expression with age via bulk RNA-seq analyses. The blue line indicates linear regression. The gray shadow covers the 95% confidence interval. Spearman correlation coefficient ρ and p values are indicated. A p value < 0.05 was considered statistically significant. n = 34 donors.

h, i. The DEGs were classified into distinct clusters by Mfuzz based on the temporal trajectories of their expression. The DEG clusters whose expression peaked in children (90 genes, h) and the elderly (144 genes, i) are shown.

j. Heatmap showing the functional enrichment analyses of the genes whose expression peaked in children in (h).

k. Heatmap showing the functional enrichment analyses of the genes whose expression peaked in the elderly in (i).

DEG, differentially-expressed genes; M-DEGs, monotonically differentially expressed genes.

Extended Data Fig. 5 | Expression dynamics of differentially expressed immune checkpoints within immune cell subsets among the different life stages.

a. Heatmaps showing the row-scaled expression of immune checkpoint ligands in the indicated cell subsets across the eight life stages. Heatmap colors from blue to red represent the relative expression level of genes. Color scale: red, high expression; blue, low expression.

b. Heatmaps showing the row-scaled expression of immune checkpoint receptors in the indicated cell subsets across the eight life stages. Heatmap colors from blue to red represent the relative expression level of genes. Color scale: red, high expression; blue, low expression.

Extended Data Fig. 6 | Global profile of T-cell clones across the entire lifespan.

a. Bar plot showing the percentage of TCR detection in each αβ T-cell subset and in total αβ T cells.

b. Bar plot showing the frequency of CDR3 amino acid length from TCR α-chains (top panel) and TCR β-chains (bottom panel) across 13 age groups.

c. Box plots showing the clonal expansion levels of each αβ T cell subset quantified by STARTRAC-expa metric in αβ T-cell subsets across 13 age groups.

CDR3, complementarity-determining region-3; TCR, T-cell receptor.

Extended Data Fig. 7 | Features of CD8_MAIT_SLC4A10 cells.

a. Scatter plots showing TCR repertoire diversity in CD8_MAIT_SLC4A10 cells across eight life stages. The diversity of TCR repertoires was measured by Inverse Simpson Index. The black line indicates LOESS regression.

b. Representative gene interaction network of upregulated genes in (Fig. 5k) using the STRING database (https://string-db.org/). Edge thickness represents the strength of data support from the STRING database.

c. Box plots showing the normalized protein expression of CD161 (encoded by KLRB1) in CD8_MAIT cells in children (1–6 y), adolescents (12–18 y), adults (30–60 y), and the elderly (70–90 y) based on the CyTOF analyses.

CyTOF, cytometry by time-of-flight; DEGs, differentially expressed genes; MAIT, mucosal-associated invariant T.

Extended Data Fig. 8 | Transcriptional regulatory characteristics of cell differentiation in the B-cell subsets.

a. FACS gating scheme for identification of the newly discovered cytotoxic B-cell subset (B_BCR+GNLY+ cells) by flow cytometry with antibodies specific for CD19, IgM, and GNLY.

b. Heatmap showing the area under the curve values of TF regulons in each B-cell subset. Numbers between brackets represent the regulon counts for respective TFs. The value for each TF regulon is a row-scaled z-score. Red, high expression; blue, low expression.

c. Rank for TF regulons in each B-cell subset by the regulon specificity score. The top 10 TF regulons are listed and marked in blue.

d. Heatmap showing the row-scaled expression of the top 10 TFs in each B-cell subset shown in (c).

FACS, representative fluorescence-activated cell sorting; TF, transcription factor.

Extended Data Fig. 9 | Functional difference between B_BCR+GNLY+ cells and other peripheral immune cells with cytotoxic functions.

a. Upset plots showing the number of shared upregulated DEGs in B_BCR+GNLY+ cells compared with that in seven cytotoxic T-cell subsets including CD4_TEM_GNLY, CD8_TCM_HAVCR2, CD8_TEM_ZNF683, CD8_TEM_CMC1, CD8_TEM_GNLY, CD8_MAIT_SLC4A10, and γδ T cells.

b. Upset plots showing the number of shared upregulated DEGs in B_BCR+GNLY+ cells compared with that in three cytotoxic NK-cell subsets including NK_CD56low, NK_CD56high, and NK_Proliferating cells.

c. Dot plots showing the representative enriched biological processes of upregulated DEGs in B_BCR+GNLY+ cells compared with those in seven cytotoxic T-cell subsets. Dot colors represent the significance of biological processes, and dot sizes indicate gene ratio annotated to each biological process.

d. Dot plots showing the representative enriched biological processes of upregulated DEGs in B_BCR+GNLY+ cells compared with those in three cytotoxic NK-cell subsets. Dot colors represent the significance of biological processes, and dot sizes indicate gene ratio annotated to each biological process.

e. Box plots showing the expression score of cytotoxicity in B_BCR+GNLY+ cells and seven cytotoxic T-cell and three cytotoxic NK-cell subsets. Box plots show the median, 25th and 75th percentiles, and whiskers extending to maximum and minimum data points. The dashed line indicates the median score of cytotoxicity in B_BCR+GNLY+ cells. P values were determined by the two-sided unpaired Wilcoxon test. A p value < 0.05 was considered statistically significant.

DEGs, differentially expressed genes.

Extended Data Fig. 10 | Features of BCR across lifespan

a. Bar plot showing the percentage of BCR detection in all peripheral immune cell subsets.

b. UMAP plots showing selected shared BCR clones among different B-cell subsets. The BCR clones are indicated in different colors.

c. Scatter plot showing the clonal expansion levels of total B cells (top left) and five B-cell subsets across eight age stages quantified by STARTRAC-expa metric. The black line indicates LOESS regression. The correlation coefficient ρ and p values are indicated. A p value < 0.05 was considered statistically significant.

d. Box plots showing the CDR3 amino acid length from BCR heavy chains (left panel) and BCR light chains (right panel) across eight life stages. After controlling for donor, Kruskal-Wallis test reveals no significant differences in CDR3 amino acid length across eight life stages (p > 0.05).

BCR, B-cell receptor; CDR3, complementarity-determining region-3; UMAP, uniform manifold approximation and projection.

Extended Data Fig. 11 | Assessment of the single-cell immune age predictive model.

a. Dot plot showing the predictive performance of the cell type-specific immune age prediction model based on our lifecycle-wide single-cell data. All cell types are shown.

b. Scatter plots showing the correlation of single-cell immune age (siAge) with chronological age (cAge) in the validation cohorts of patients with disturbed immune function: nasopharyngeal carcinoma (NPC) cohort (left panel), Kawasaki diseases (KD) cohort (middle panel), and systemic lupus erythematosus (SLE) cohort (right panel). The blue line indicates linear regression. The gray shadow covers the 95% confidence interval. Pearson correlation coefficient (PCC) and p values are indicated. A p value < 0.05 was considered statistically significant.

siAge: Model to Predict the Immune-Age

Single Cell Online Explorer