Human autoimmunity at single cell resolution in aplastic anemia before and after effective immunotherapy

Abstract

Severe immune aplastic anemia is a fatal disease due to the destruction of marrow hematopoietic cells by cytotoxic lymphocytes, serving as a paradigm for marrow failure syndromes and autoimmune diseases. To better understand its pathophysiology, we apply advanced single cell methodologies, including mass cytometry, single-cell RNA, and TCR/BCR sequencing, to patient samples from a clinical trial of immunosuppression and growth factor stimulation. We observe opposing changes in the abundance of myeloid cells and T cells, with T cell clonal expansion dominated by effector memory cells. Therapy reduces and suppresses cytotoxic T cells, but new T cell clones emerge hindering robust hematopoietic recovery. Enhanced cell-cell interactions including between hematopoietic cells and immune cells, in particular evolving IFNG and IFNGR, are noted in patients and are suppressed post-therapy. Hematologic recovery occurs with increases in the progenitor rather than stem cells. Genetic predispositions linked to immune activation genes enhances cytotoxic T cell activity and crosstalk with target cells.

Fig. 1. Resolving hematopoietic and immune cell complexity using CyTOF and scRNA-seq.

a Experimental workflow. BMMNC samples from patients and healthy donors (HD) were analyzed by CyTOF and flow cytometry to profile hematopoietic and immune cells and subpopulations. BMMNCs and Lineage⁻CD34⁺ cells were subjected to scRNA-seq, followed by data analysis including single-cell transcriptome profiling (gene expression and cell–cell interaction) and scTCR/BCR profiling. Created in BioRender. Wu, Z. (2025) https://BioRender.com/t6tfzun. b Characterization of the patient cohort. c UMAP projection of BMMNCs at all time points from patients and healthy donors in CyTOF. d Log fold changes and log P values (with the two-sided unpaired Mann–Whitney test) of cell abundance between pre-treatment samples of SAA patients and controls (top) and post-treatment samples (bottom). Red dots represent lymphoid cells, and blue dots represent myeloid cells. e UMAP embedding of CyTOF data as shown in Fig. 1c, colored by differential abundance. Red indicates clusters with patients’ cells (pre-treatment samples) constituting >95% of cells, and blue clusters with healthy donors’ cells >95%. f Graphic of Nhoods identified by Milo. Nodes are Nhoods, colored by log2 fold-changes (log2FC) between pre-treatment samples of SAA patients (n = 20) and healthy donors (n = 4). Nondifferential abundance Nhoods (FDR ≥ 0.1) are colored white, and sizes correspond to the number of cells in a Nhood. Graph edges depict the number of cells shared between adjacent Nhoods. g beeswarm plot showing adjusted log2FC of cell populations abundance between pre-treatment samples of SAA patients and healthy donors in Nhoods. h A UMAP plot of single-cell gene expression of BMMNCs samples at different time points of patients (n = 45) and healthy donors (n = 4). Cells are colored by type (HSPC, ProB as B cell progenitors, erythroblast, neutrophil lineage, monocytic lineage, CD8⁺ T cell, CD4⁺ T cell, NK cell, B cell and plasma cell, dendritic cell, stromal cell, and platelet). i Frequency (% in BMMNCs) of CD8⁺ T, CD4⁺ T, NK, erythroblast, neutrophil, and monocyte in pre- (n = 20), post-treatment (n = 14) samples of SAA patients and healthy donors (n = 4). P values with the two-sided unpaired and paired Mann–Whitney test are shown.

Fig. 2. Loss of early-stage and recovery of late-stage HSPCs after treatment.

a Representative plots of phenotypes of HSPCs in healthy donors (HD) and patients by flow cytometry. b Proportions of progenitor populations were compared between pre- (n = 20) and post-treatment samples (n = 16) and HD (n = 4). P values with the two-sided unpaired and paired Mann–Whitney test. c A UMAP plot of single-cell gene expression in HSPCs of all patients and HD, colored by cell types as HSC, MEP, GMP, Tprog, and Bprog. d A heatmap of representative differentially expressed genes (grouped by functional pathways) of HSPCs between patients (n = 20) and HD (n = 4) on the top panel, between post- (n = 16) and pre-treatment samples (n = 20) at the bottom. Values are presented as log2FC. e A principal component analysis plot of samples based on cell cluster abundance measured by CyTOF. Pre- and post-treatment samples and samples from HD were clustered separately. f A plot showing fraction of variance (gene expression of BMMNCs) explained by diseases (PT/HD), disease severity (SAA/VSAA), and treatment (Post/Pre-immunosuppression). g Dynamic changes of lineage priming of HSCs to LymP, GMP, and MEP, along with pseudotime differentiation. X-axis, pseudotime ordering from HSCs to lineage-restricted progenitors estimated by Palantir. Y-axis, Log(lineage signature gene expression) in pre-treatment (dark orange), post-treatment samples (orange), and HD (gray). Data are presented as mean with 1.96SE. P values with the Analysis of Covariance (ANCOVA). h Dynamic changes of expression levels of transcription factors along pseudotime differentiation. Y-axis, expression of transcription factors in pre-treatment (dark orange), post-treatment samples (orange), and healthy donors (gray). Data are presented as mean with 1.96SE. P values with ANCOVA. i Dynamic changes of lineage primed HSC fraction in pre-treatment (dark orange), post-treatment samples (orange), and HD (gray) along differentiation. Y-axis, Log(cell numbers of HSCs). For each lineage, cell fractions of early (pseudotime 1–4 on x-axis) and late stages (pseudotime 5–10) were compared: a ratio in pre-treatment/HD, and a ratio in post-treatment/pre-treatment. P values with the two-sided unpaired t-test: HSC to GMP, 0.016 and 0.007; HSC to MEP, 0.929 and 0.96; HSC to LymP, 0.004 and 0.029.

Fig. 3. CD8⁺ T cells expanded and exhibited cytotoxicity in SAA.

a In a UMAP only showing CD8⁺ T cells, subtypes including naïve, central memory (CM), and effector memory (EM) were indicated. Pseudotime based on Slingshot was calculated and colored accordingly. b RNA velocity analysis of CD8⁺ T cell phenotypes independently confirmed the differentiation trajectory. c Dot plots showing frequency of naïve, CM, and EM subpopulations of CD8⁺ T cells in pre- (n = 20), post-treatment (n = 18) samples of SAA patients and healthy donors (n = 4). P values with the two-sided unpaired and paired Mann–Whitney test are shown. d Plots showing expression dynamics of marker or functional genes of CD8⁺ T cells in pre- and post-treatment samples in SAA patients and healthy donors, along pseudotime. e Plots of frequency of effector memory T cells and cytotoxicity gene expression of CD8⁺ T cells in pre- and post-treatment samples of SAA patients and healthy donors.

Fig. 4. CD8⁺ T cells and CD4⁺ T cell subpopulations expanded and activated in SAA.

a UMAP plots showing overlay of CD8⁺ T cell density in healthy controls (n = 4), SAA patients pre- (n = 20) and post-treatment (n = 18), as well as the individual conditions. Dot plots showing abundance (% in CD45⁺ BMMNCs) of CD45RA⁺CD8⁺CD127⁺CD197⁺ naïve T cells (cluster 4), CD45RA⁺CD8⁺CD57⁺ effector T cells (cluster 5), and CD45RA⁺CD197⁺CD184⁺CD127⁺ central memory T cells (CM, cluster 15) in pre- and post-treatment samples of SAA patients and healthy donors identified by Binary Clust. P values with the two-sided unpaired and paired Mann–Whitney test are shown. b UMAP plots showing overlay of CD4⁺ T cell density in healthy donors (n = 4), SAA patients pre- (n = 20) and post-treatment (n = 18) as well as individual conditions. Dot plots showing abundance (% in CD45⁺ BMMNCs) of CD45RA⁺ CD4⁺CD197⁺CD38⁺ naïve T cells (Cluster 18), CD45RA⁺CD4⁺CD197⁺CD38⁺CD127⁺CD27⁺CD184⁺ central memory T cells (CM, cluster 16), CD4⁺CD25⁺CD127⁺CD45RO⁺CD150⁺FOXP3⁺CD95⁺CD194⁺CD39⁺ Treg-B (cluster 8), CD4⁺CD25⁺CD127^-CD45RA⁺Foxp3⁺CD197⁺ Treg-A (cluster 6), and nonA nonB-Treg (cluster 11) in pre- (n = 20), post-treatment (n = 18) samples of SAA patients and healthy donors (n = 4) identified by BinaryClust. P values based on the two-sided unpaired and paired Mann–Whitney tests. c Ratios of Treg-A and Treg-B abundance versus CD8⁺ effector T cells (cluster 5 and cluster 2) in pre- (n = 20), post-treatment (n = 18) samples of SAA patients and healthy donors (n = 4). P values based on the two-sided unpaired and paired Mann–Whitney tests. d In a UMAP only showing CD8⁺ T cells, CD57⁺ T cells, and expression of cytotoxicity genes were indicated. e In a UMAP only showing CD4⁺ T cells, Treg-A and Treg-B were shown. f In a UMAP only showing CD4⁺ T cells, subtypes including naïve, central memory (CM), and effector memory (EM) were indicated. A bar chart showing percentages of Treg-B and non-Treg-B CD4⁺ T cells as naïve, CM, and EM. g Gene Set Enrichment Analysis (GESA) enriched plots of differentially expressed genes for Hallmark interferon gamma response and Hallmark interferon alpha response in Treg-B compared with non-Treg-B CD4⁺ T cells in pre-treatment samples of SAA patients. GSEA is based on the one-sided Kolmogorov–Smirnov test. h A heatmap showing expression of representative differentially expressed genes grouped by their functional pathways in IFN-γ and IFN-α signaling between Treg-B and non-Treg-B CD4⁺ T cells pre- and post-treatment in SAA patients. Values are presented as log2FC.

Fig. 5. Effector memory CD8⁺ T cells were clonally expanded in SAA.

a Frequency of clonal expanded T cells in pre- (n = 20), post-treatment (n = 16) samples, and HD (n = 4). P values with the two-sided unpaired and paired Mann–Whitney test. Clone sizes were projected to a UMAP of T cells in pre-treatment samples, colored based on clone sizes. b Gini index of TCR clone sizes in pre- (n = 20), post-treatment (n = 16) samples, and HD (n = 4). P values with the two-sided unpaired and paired Mann–Whitney test. c CD4⁺ T and CD8⁺ T cells expressing the highest (top 10%) cytotoxicity score and IFN-γ signaling score are highlighted in red, and all the rest in gray. Expression of T cytotoxicity genes and TNF-α via NFκB signaling genes was compared between clonally and non-clonally expanded T cells. P values with the two-sided unpaired and paired Mann–Whitney test. Bar chart showing percentage of clonally expanded and nonexpanded T cells as naïve, CM, and EM. d Expression dynamics of gene scores in averaging increased, decreased, and stable clones after treatment. P values were generated using the two-way ANOVA with factors pre/post and increase/decrease/unchanged clones.

Fig. 6. T cell clonal expansion dynamics associated with hematopoietic recovery.

a Patients were grouped: clone sizes increased (red, n = 7) and decreased (blue, n = 9). Data are presented as mean values ± SEM. P values with the two-sided unpaired Mann–Whitney test were shown. b Clinical characteristics of patients with clone sizes increased (n = 7) and decreased (n = 9) post-treatment. P values with the two-sided unpaired Mann–Whitney test are shown. c New clone sizes in two groups of patients who had clone sizes increased (n = 7) and decreased (n = 9) were compared. P values with the two-sided unpaired Mann–Whitney test are shown. Correlation of new clone sizes with total clone size changes after treatment. P values and slopes with the Pearson correlation test are shown. d Clone sizes (frequency %) of pre-existing clones in paired samples of patients (n = 16) and HD (n = 4). P values with the two-sided unpaired and paired Mann–Whitney test are shown. e Correlation of pre-existing clone sizes with new clone sizes (left) and new clone sizes with therapeutic scores (right). P values and slopes with the Pearson correlation test are shown. f TCR sequence similarity is plotted for T cells in HD, pre-existing clones in pre-treatment samples (pre), new clones in post-treatment samples (post), and between pre-existing and new clones in post-treatment samples (post). Those who had higher similarities of pre-existing clones are highlighted in red across sample groups. P values with the two-sided unpaired Mann–Whitney test are shown. g Expression of IFN-γ signaling scores was plotted for all patients (n = 20), for patients who had clone sizes increased (n = 7) and decreased (n = 9), and HD (n = 4), respectively. P values with the two-sided unpaired and paired Mann–Whitney test are shown.

Fig. 7. Enhanced cell–cell interactions in SAA and attenuation by therapy.

a Ligand-receptor pairs among cell types in BMMNCs were estimated by CellPhoneDB. Color legends for cell types are the same as in Fig. 2a. The Thickness of lines connecting cell types indicates the total number of ligand-receptor pairs between two cell types estimated by CellPhoneDB. In general, there were more ligand-receptor interactions between cell types in BMMNCs of SAA (left) than in those of healthy donors (right). CD8T CD8⁺ T cell, CD4T CD4⁺ T cell, NK natural killer cell, CD34 CD34⁺ cell, ProB proB cell, B_Plasma B cell_Plasma cell, Ery erythroblast, Neut neutrophil, Mono monocyte, DC dendritic cell. b Heatmaps showing LogFC of cell–cell interaction scores estimated by CellPhoneDB across cell types in BMMNCs post-treatment as compared to pre-treatment. c Ligand-receptor pairs that were overrepresented in SAA patients than in healthy donors across cell types in the BMMNCs are presented, with HSPCs as sender (left) and receiver (right). Significance indicates if the ligand-receptor pair is overrepresented in pre-treatment samples compared with controls. d Cell–cell interactions of IFNG/IFNGR. Heatmaps depicting the interactions of IFNG and IFNGR1/2 between CD8⁺ T cells (left, as sender and right, as receiver). A sum of IFNG/IFNGR interactions in CD8⁺ T cells with all other cell types in the BM is shown at the bottom for healthy donors, pre-, and post-treatment samples. e A heatmap depicting interactions of IFNG and IFNGR1/2 between various cell types in the BM (as sender) and HSPC (as receiver). Sum of IFNG/IFNGR interactions in HSPCs with all other cell types in the BM is shown at the bottom for healthy donors, pre-, and post-treatment samples. f Heatmaps showing ligand-receptor pairs between HSPCs and pre-existing T cell clones in pre-treatment samples (left), and between HSPCs and residual pre-existing T cell clones, and HSPCs and new T cell clones (right). g Heatmaps showing numbers of ligand-receptor pairs between CD8⁺ T cells and early-stage HSCs, late-stage HSCs, HSCs, and progenitors.

Fig. 8. Genetic risk and cell-type-specific disease-related variants contribute to disease phenotype.

a A Manhattan plot of GWAS analysis highlighting risk genetic loci for AA, based on GCST9001879 data, in which P values were calculated by the generalized linear mixed model with adjusting for age, sex, and top 20 principal components. b A Circus plot showing top risk genes of AA of MAGMA analysis. The outer ring demonstrates 22 autosomal chromosomes. In the inner ring, a circular symbol represents a specific risk gene. c GESA enriched plots of risk genes identified in GWAS analysis were enriched in Hallmark interferon gamma response. GSEA based on the one-sided Kolmogorov–Smirnov test. d UMAP embedding BMMNCs from healthy donors (left) and SAA patients (right), colored by scDRS disease scores calculated from GWAS summary statistics. e A dot plot showing results of combination of scRNA-seq data and GWAS summary statistics on AA based on RolyPoly analysis among patients. Y-axis shows mean negative log-transformation P value (-Log2P), x-axis shows cell types. P values (one-sided) were calculated by estimating the association of each cell type with a phenotype. f A heatmap showing percentages of cells expressing selected risk genes in multiple cell types. g Same UMAP projection of BMMNCs as in Fig. 2a. Expression of S1PR5 and IRF5 is plotted, showing cells (in red) with top 10% expression levels of S1PR5 and IRF5. h Frequency of S1PR5⁺CD8⁺ T cells and S1PR5⁺NK cells were compared between pre-treatment (n = 20) and post-treatment samples (n = 16) of SAA patients and healthy donors (n = 4). P values with the two-sided unpaired and paired Mann–Whitney test are shown. i Immune activation scores in S1PR5⁺ and S1PR5⁻ CD8⁺ T cells, S1PR5⁺ and S1PR5⁻ NK cells, and IRF5⁺ and IRF5⁻ myeloid cells. P values with the two-sided unpaired and paired Mann–Whitney test are shown. j Top downregulated gene pathways in NFκB⁺ compared with NFκB⁻ myeloid cells. k Cell interaction scores calculated by CellPhoneDB in S1PR5⁺ and S1PR5⁻CD8⁺ T cells, NK cells, IRF5⁺ and IRF5⁻ myeloid cells, and NFκB⁺ and NFκB⁻ myeloid cells.