Demographic and clinical characteristics of individuals with AA

Table 2 provides an extended list of the clinical information of the individuals included in the analysis.

Table 2 Clinical characteristic of healthy donors and patients with AAAlters of immune profile of AA patients

To broadly profile the individual components of the immune response against aplastic anemia and to assess the general landscape of immune responses and their perturbation, we performed extensive immunophenotyping to characterize the frequencies of circulating immune subsets in NSAA, SAA and VSAA compared to HDs, we collected longitudinal PBMCs and serum samples from adults enrolled in Guangzhou First People’s Hospital and performed flow cytometry analysis. We employed an antibody panel targeting 27 cellular markers to identify cell types and assess immune activation, function, and proliferation. Due to the proclivity of FlowJo software to encounter crashes when processing voluminous datasets for t-SNE downscaling, we judiciously selected 9 samples per subgroup for this analysis: HDs (n = 9), NSAA (n = 9), SAA (n = 9), and VSAA (n = 9). This selective methodology is an established standard in flow cytometry, intended to mitigate potential biases and anomalies that could compromise the integrity of data interpretation. After preprocessing, cells were clustered at three resolutions: five main types (B cell, natural killer (NK) cell, myeloid, monocyte, and T cell), seven cell subtypes. Lineage-specific protein expression guided cell clustering to minimize bias from activation, function, and proliferation proteins. T-distributed stochastic neighbor embedding (t-SNE) analysis revealed differential abundance of immune cell phenotypes across the entire cohort and various disease severity categories, irrespective of cluster labels (Fig. 1A-B). We employed FlowJo software to manually gate flow cytometry data (Supplementary Fig. 4), defining distinct cell populations. This formed the basis for subsequent analyses. To visualize differences in cell population abundances between HDs and AA patients, we performed t-SNE dissect the difference between HDs, NSAA, SAA and VSAA (Fig. 1C). We observed a decline in the proportion of LDNs populations in VSAA compared to HDs (p = 0.0126 for LDNs, Fig. 1D). The LDNs frequency also differed between NSAA and VSAA (p = 0.0289). The abundance of eosinophils was significantly greater in NSAA and SAA compared to HDs. In contrast, the abundance of monocytes was significantly lower in SAA and VSAA compared to NSAA (p = 0.0012 for SAA, p = 0.0332 for VSAA). Furthermore, we did not find significant differences in the frequencies of natural killer T cells (NKT) and NK cells between HDs and NSAA, SAA or VSAA (Fig. 1D). The application of Principal Component Analysis (PCA) and dimensionality reduction to the entire dataset of immune cell data from AA and HDs demonstrated a clear separation between the two groups, signifying substantial differences (Fig. 1E). This finding confirms that our analytical framework effectively captured significant variations in the immune state. Collectively, these results demonstrate the efficacy of our experimental and analytical approaches in identifying variations in the abundance of immune cell subtypes that are associated with different levels of disease severity in AA.

Fig. 1

Alters of immune profile of AA patients. A Experimental outline showing peripheral blood mononuclear cells (PBMC) collection and mass cytometry analysis. B T-distributed stochastic neighbor embedding (t-SNE) plot of major immune cells in PBMCs. Cells are colored based on cell types; heat map of major immune cells in PBMCs, clustered by their relative expression of the markers. C T-SNE projections of major immune cells in PBMCs: the healthy Donors (HDs), non-severe aplastic anemia (NSAA) groups, severe aplastic anemia (SAA)groups and very severe aplastic anemia (VSAA) groups, respectively. D Percentage of major immune cells in PBMC from the HDs (n = 9), NSAA group (n = 9), SAA group (n = 9), VSAA group (n = 9). In all plots, mean ± SD is shown, and p > 0.05 was no statistically significant difference. E Utilizing PCA to reduce the dimensionality of the data, analyze the differences between HDs and AA

NK cells exhibit a less functional CD56 Dim molecular signature in AA patients

In the following analysis, we concentrated on distinct immune cell classes and utilized dimensionality reduction methods. Figure 2A illustrates the downscaled map obtained through clustering for all NK cells. These NK cells can be further classified into distinct subgroups based on the marker heat map, which identifies the following categories: NKT, NKBri (immature NK cells), NKDim (mature NK cells), and NKUC (unclassical NK cells) (Ming et al. 2020). To elucidate variations in cell population proportions between HDs and patients with AA, t-SNE analysis was employed to identify distinct patterns among HDs, NSAA, SAA, and VSAA (Fig. 2B). This visualization technique facilitates the exploration of high-dimensional data, allowing us to uncover and interpret the differences in cell populations across these diverse conditions. The figure indicates that the cellular clustering patterns of various subtypes vary within each group. Various subtypes of NK cells were classified using manual Flow Cytometry gating (Fig. 2C). Compared to the HDs, we observed that NKDim cells exhibit a significant decrease in both the SAA and VSAA groups, while NKUC cells show a noticeable increase in the NSAA and VSAA groups. NKBri cells do not show significant changes across all groups (Fig. 2D). Research has demonstrated that NKBri cells serve as precursors for NKDim cells, with the latter exhibiting shorter telomeres compared to the former. NKBri cells primarily function in secreting large amounts of cytokines, while NKDim cells play a role in exerting cytotoxic effects (Poli et al. 2009). The majority of CD56+ NK cells express CD27. During the transition from CD56Bri to CD56Dim cytotoxic effectors, there is a reduction in CD27 expression (Vossen et al. 2008). Consequently, a statistical analysis of CD56−CD27−cells within NK cells was conducted. Compared to the HDs group, a significant decrease was observed in the AA group, consistent with the changes observed in NKDim cells in AA (Fig. 2E). This suggests a significant reduction in less functional CD56Dim NK cells in AA. In conjunction with clinical data, we also observed a negative correlation between the proportion of CD56Dim cells in NK cells and C-reactive protein (CRP) in AA, aligning with the afore-mentioned results (Fig. 2F). Subsequently, we examined the expression levels of relevant functional factors (NKG2A, NKG2C, NKG2D, NKp30, NKp44, and NKp46) in NKDim cells. A decrease in expression levels was found in the AA group compared to the healthy group (Fig. 2G), consistent with previous studies (Alter et al. 2004). Additionally, there was a significant decrease in the expression levels of the cytotoxicity factor CD107a in NK cells and each subtype compared to the HDs group (Fig. 2H). In summary, there is a significant reduction in less functional CD56Dim NK cells in AA, indicating functional impairment.

Fig. 2

NK cells exhibit a less functional CD56Dim molecular signature in AA patients. A t-SNE plot of NK cells of PBMCs from all samples; heat map of major cell subtypes in NK cells, clustered by their relative expression of the markers. B t-SNE plot of major NK cell subsets. Cells are colored based on cell types. C The flow cytometry gating strategy for NK cell phenotyping. D Percentage of major NK cell subsets in NK cells from the HDs (n = 9), NSAA group (n = 9), SAA group (n = 9), VSAA group (n = 9). E The flow cytometry gating strategy for CD56Dim NK cells, as well as the statistical plots for CD56DimCD27−NK cells and CD56BriCD27−NK cells. F The correlation between the ratio of CD56Dim NK cells to CD56+ NK cells and C-reactive protein (CRP). p < 0.05, it’s statistically significant. G Statistical plots of the expression levels of relevant functional factors (NKG2A, NKG2C, NKG2D, NKp30, NKp44, and NKp46) in NKDim cells in the AA and HDs groups. p < 0.05 are considered statistically significant and labeled accordingly. H Statistical plots of the expression levels of the cytotoxic factor CD107a in NK cells and each subtype in the HDs and AA groups. In all plots, mean ± SD is shown, and p < 0.05 are considered statistically significant and labeled accordingly

CD56 +(NK-liked) monocytes increase in AA patients and possess NK peculiarity

We employed t-SNE plots to visualize the dimensionality reduction outcomes of various mononuclear cell subtypes, including classical monocytes (c-monocytes), intermediate monocytes (inter-monocytes), non-classical monocytes (non-c monocytes), and CD56+ monocytes. These subgroups were delineated using flow cytometry marker heatmaps (Fig. 3A). Next, we categorized the samples into four groups: HDs, NSAA, SAA, and VSAA. Following dimensionality reduction analysis, we observed differences in the clustering patterns of each mononuclear cell subtype among these groups (Fig. 3B). Figure 3C depicts the flow cytometry profiles of different mononuclear cell subtypes. The proportions of c-monocytes, inter-monocytes, and non-c monocytes among mononuclear cells did not significantly differ across HDs, NSAA, SAA, and VSAA groups. Nevertheless, CD56+ monocytes showed a significant increase in NSAA and VSAA compared to HDs (Fig. 3C). Using manual gating in flow cytometry, we detected a substantial increase in CD56+ monocytes in AA, consistent with findings in rheumatic diseases and post-COVID-19 infection, suggesting a contributing role (Kuri-Cervantes et al. 2020). To further analyze the expanded CD56+ monocyte population in t-SNE plots, we utilized manual gating in flow cytometry (Fig. 3D). Additionally, ROC curve analysis indicated that the alterations in CD56+ monocytes are specific to AA (Fig. 3E, p < 0.0001). Furthermore, within CD14+ monocytes in AA patients exhibited a significant increase in CD56 expression (Fig. 3F). Additionally, in healthy individuals, the CD56+CD14+ cell subset displays NK cell-like attributes, characterized by GZMB, Perforin, and T-bet expression, indicating cytotoxic properties. Conversely, these functions are notably impaired in AA patients (Fig. 3G). In summary, the innate immune cells in AA patients, including NK cells and NK-like monocytes, both exhibit defects in natural killing and phagocytic functions.

Fig. 3

CD56+ (NK-liked) monocytes increase in AA patients and possess NK peculiarity. A t-SNE plot of monocytes of PBMCs from all samples; heat map of major cell subtypes in monocytes, clustered by their relative expression of the markers. B t-SNE plot of major monocyte subsets. Cells are colored based on cell types. C percentage of major monocyte subsets in monocytes from the HDs (n = 9), NSAA group (n = 9), SAA group (n = 9), VSAA group (n = 9). D The proportion of CD56+ monocytes within the monocyte population in the HDs and AA groups, along with statistical plots. E ROC curve analysis of CD56+monocyte in AA. F The difference in CD56 expression within CD14+ monocytes between AA and HDs. G The expression differences of GZMB, Perforin, and T-bet on NK cells between AA and HDs. In all plots, mean ± SD is shown, and p < 0.05 are considered statistically significant and labeled accordingly

MDSCs decrease in AA and recover post-treatment

We evaluated the proportions of MDSCs in the peripheral blood (PB) of AA patients with different subtypes (NSAA, SAA, and VSAA) and healthy donors (HDs). Utilizing previously published methods (He et al. 2018), we marked MDSCs with CD15, CD14, HLA-DR, and CD11b surface markers and visualized the flow cytometry data using t-SNE analysis (Fig. 4A-B). For further subtype analysis, MDSCs were identified with the specific markers HLA-DR−CD33+CD11b+ and classified into CD15+ polymorphonuclear MDSCs (PMN-MDSCs), CD14+ monocytic MDSCs (M-MDSCs), and CD15−CD14− early MDSCs (e-MDSCs) (Fig. 4C). Previous research has pointed to a significant decrease in MDSC proportions in AA, linked to a noticeable reduction in immunosuppressive capabilities (Dong et al. 2022). Our study similarly observed a gradual decrease in MDSC proportions in AA patients' PB compared to HDs, with no significant differences detected among the subtypes (Fig. 4D). We further analyzed samples from AA patients who received anti-thymocyte globulin (ATG) immunosuppressive therapy, dividing them into recovery (Re) and non-recovery (Non-Re) groups (Table 3). Results showed a significantly higher frequency of MDSCs in the recovery group compared to the non-recovery group, with an overall increase post-treatment (Fig. 4E). These findings suggest MDSCs as a potential marker for treatment evaluation. Using the same blood samples for both complete blood count (CBC) and flow cytometry staining enabled us to evaluate the correlation between MDSCs and various hematological indicators. Consequently, the results demonstrated a significant association between MDSC levels and key hematological parameters, specifically platelet count (PLT) and white blood cell count (WBC) (Supplementary Fig. 5). To uncover the molecular mechanisms behind the reduced suppressive function of MDSCs in AA, we performed Smart-seq analysis on MDSCs isolated from PBMCs. PCA analysis revealed distinct gene expression profiles of MDSCs between HDs and AA patients (Fig. 4F). Volcano plots and KEGG analysis highlighted a marked upregulation of the one-carbon synthesis pathway in MDSCs from AA patients (Fig. 4G). Additionally, GSEA analysis demonstrated diminished inhibitory function in AA patient MDSCs, associated with pathways such as JAK/STAT signaling, endoplasmic reticulum stress, and graft-versus-host disease (Fig. 4H-I).

Fig. 4

MDSCs decrease in AA and recover post-treatment. A t-SNE plot of MDSCs from all samples; heat map of major cell subtypes in MDSCs, clustered by their relative expression of the markers. B t-SNE plot of major MDSC subsets. Cells are colored based on cell types. C The more specific marker gating strategies for detailed analysis of MDSC subtypes, statistical chart plotting, mass spectrometry analysis, and clinical correlation analysis. D Percentage of major MDSC subsets in MDSCs from the HDs (n = 9), NSAA group (n = 9), SAA group (n = 9), VSAA group (n = 9). E Percentage of major MDSC subsets in CD45+ cells from the HDs, Acquired Aplastic Anemia (AA) group, recovery (Re)group, non-recovery (Non-Re) group. F Utilizing PCA to reduce the dimensionality of the data, analyze the differences of MDSCs between HDs and AA. G The volcano plots and KEGG analysis revealed alterations in the one-carbon synthesis pathway in MDSCs derived from AA. H GSEA analysis revealed the functionality of MDSCs in AA patients. I The enrichment plot illustrates the pathways associated with the functional changes of MDSCs in AA patients. In all plots, mean ± SD is shown, and p < 0.05 are considered statistically significant and labeled accordingly

Table 3 Hematologic recovery data tableMDSCs effectively discriminate between acquired aplastic anemia and congenital aplastic anemia

CAA primarily results from genetic mutations, with the diagnosis of congenital bone marrow failure relying mainly on NGS sequencing. In AAA, MDSCs exert immunosuppressive effects, leading to a significant decrease in their components (Dong et al. 2022). To further investigate the differences between CAA and AAA, we collected data from patients diagnosed with either CAA or AAA (Table 4). We evaluated the percentage of HLA-DR−CD11b+CD33+ MDSCs in AAA and CAA, revealing a marginal decrease in MDSCs in AAA compared to the HDs. Conversely, CAA demonstrated a pronounced increase in MDSCs (Fig. 5A). Next, the frequencies of the three MDSCs subtypes (PMN-MDSCs, M-MDSCs, and e-MDSCs) were further examined in CD45+ cells in HDs, AAA and CAA. It was observed that all three MDSCs subtypes exhibited significant differences in frequencies between AAA and CAA, proving to be efficient markers for distinguishing between the two conditions (Fig. 5B). Through ROC curve analysis, we illustrated the distinctions in MDSC populations among the HDs, AAA, and CAA groups. Notably, significant differences were evident between the control group and both the AAA and CAA groups, as well as between AAA and CAA (Fig. 5C). Subsequently, as potentially effective indicators for clinical application, we successfully translated and implemented these findings in a clinical context. In CAA patients, there exists a primary association between telomere length and the extent of chromosomal breakage (Gramatges and Bertuch 2013). Our observations indicate that CAA patients with shorter telomeres exhibit elevated levels of MDSCs (Fig. 5D), with a negative correlation observed between telomere length and the proportions of MDSC subtypes (Fig. 5E). Furthermore, patients with higher degrees of chromosomal breakage in CAA demonstrate increased proportions of MDSCs (Fig. 5F), revealing a significant positive correlation between the degree of chromosomal breakage and MDSC proportions (Fig. 5G). Collectively, our findings highlight the distinguish roles of MDSCs in distinguishing between acquired and congenital aplastic anemia.

Table 4 Clinical characteristic of healthy donors and patients with AAA and CAAFig. 5

MDSCs effectively discriminate between AAA and CAA. A Utilizing flow cytometry to discern variances in MDSCs between the healthy control group, AAA, and CAA. B percentage of major MDSC subsets in CD45+ cells from the control group, AAA group, CAA group. C The ROC curves demonstrate the correlations of MDSCs between the control group, AAA and CAA pairwise. D The respective proportions of MDSCs with shorter telomeres and MDSCs with longer telomeres. E The correlation between the subtypes of MDSCs and telomere length. FThe proportions of MDSCs with chromosomal breakage over 10% and under 10% respectively, among the total MDSC population. G The correlation between the subtypes of MDSCs and chromosomal breakage. In all plots, mean ± SD is shown, and p < 0.05 are considered statistically significant and labeled accordingly

AA patients exhibit elevated activated cytotoxic T cells and reduced regulatory T cell

We conducted dimensionality reduction on flow cytometry data and visualized the results using t-SNE. This approach aimed to examine the differential expression levels of T-cell subtypes between HDs and individuals with AA (Fig. 6A-B). In our study, we found no significant differences in various subsets of CD4+ and CD8+ T cells-including naive T cells (Tnaive), central memory T cells (Tcm), effector memory T cells (Tem), and effector memory cells re-expressing CD45RA (Temra) between AA patients and HDs (Fig. 6C-D). Previous research has highlighted the crucial role of T cells in AA (Meraviglia et al. 2019). Compared to HDs, AA patients exhibit elevated levels of cytotoxic CD8+ T cells in both bone marrow and PB, coupled with a decrease in Treg levels (Wang and Liu 2019). Furthermore, we observed a reduced proportion of Tregs within the CD4+ T cell subset in the AA compared to HDs. Additionally, we conducted supplementary assessments of Treg activation status by examining the surface expression of CD45RA, revealing a lower abundance of activated Tregs in the peripheral blood of AA patients. This confirms both a reduction in Treg levels and functional impairment in AA (Fig. 6E). Moreover, our investigations found no significant differences in the proportion of activated CD4+ T cells between HDs and AA. However, in comparison to HDs, there was a higher proportion of activated CD8+ T (CD38+HLA-DR+) cells in VSAA (Fig. 6F). Prior studies have suggested that activated CD8+ T cells induce inflammation (Kuri-Cervantes et al. 2020). Additionally, we employed t-SNE for dimensionality reduction and visualization to explore variations in B cells and their subsets in AA (Supplementary Fig. 6A-B). No significant differences were observed in the three B cell subtypes among NSAA, SAA and VSAA (Supplementary Fig. 6C), consistent with previous research (Kordasti et al. 2012). In summary, our findings suggest alterations in T cell subsets in AA, characterized by an increase in cytotoxic T cells, a concurrent decrease in Treg cells, and functional dysregulation.

Fig. 6

AA exhibit elevated activated cytotoxic T cells and reduced regulatory T cells. A t-SNE plot of T cells from all samples; heat map of major cell subtypes in T cells, clustered by their relative expression of the markers. B t-SNE plot of major T cell subsets. Cells are colored based on cell types. C Percentage of major CD4+ T cell subsets in MDSCs from the HDs (n = 9), NSAA group (n = 9), SAA group (n = 9), VSAA group (n = 9). D Percentage of major CD8+ T cell subsets in MDSCs from the HDs (n = 9), NSAA group (n = 9), SAA group (n = 9), VSAA group (n = 9). E The level and functional alteration of regulatory T cells in AA. F The changes in the proportion of activated T cells within the T cell population across the HDs (n = 9), NSAA group (n = 9), SAA group (n = 9), VSAA group (n = 9). In all plots, mean ± SD is shown, and p < 0.05 are considered statistically significant and labeled accordingly

Immune cell correlation identification during AA

Based on the comprehensive immune mapping, we utilized three statistical metrics—log fold change (Log FC), area under the curve (AUC), and false discovery rate (FDR)—to analyze the differences in immune cells HDs and AA. As shown in Fig. 7A, memory B cells, eosinophils, and NKBri cells were significantly elevated in AA (Log FC > 0.25, FDR < 0.05, AUC > 0.7), whereas NKDim cells and non-classical monocytes were significantly reduced (Log FC < -0.25, FDR < 0.05, AUC > 0.7). The results from the volcano plot were consistent with those from the comprehensive immune mapping, further validating the scientificity.

Fig. 7

Immune cell correlation and novel immune cells identification during AA. A Plot of correlation between HDs and AA according to Log FC and -log10 (FDR); Plot of correlation between HDs and AA according to Log FC and AUC. B The correlation heatmap of various key immune cells in AA. C Displaying the correlation among various immune cells using a Chord Diagram

By analyzing multiple immune cell types within the same samples, we were able to more scientifically assess the relationships among these cells, providing new insights for subsequent research. We employed Spearman's correlation test to evaluate the correlations among major immune cells and their subtypes (Fig. 7B), and combined with the chord diagram (Fig. 7C). This approach revealed significantly correlated cell pairs (p < 0.05). Specifically, NKBri cells showed significant positive correlations with monocytes, NKT cells, various CD8+ T cell subsets, B cells, and c-monocytes. NKDim cells were positively correlated with B cells, inter-monocytes, NKT cells, activated CD4+ T cells, LDNs, and non-c monocytes. C-monocytes exhibited significant positive correlations with B cells, inter monocytes, and other monocyte subsets, while naive B cells were positively correlated with memory B cells (p < 0.05). Additionally, NKBri cells were negatively correlated with CD4+ Tcm, and NKDim cells showed negative correlations with CD4+ Temra, NKUC, CD4+ Tem, and naive B cells. Eosinophils were negatively correlated with regulatory T cells, whereas classical monocytes exhibited significant negative correlations with various T cell subsets (p < 0.05). These findings suggest that by leveraging Comprehensive mapping of immune, we can delve deeper into the significant alterations in immune cells in AA, exploring their mechanisms and interactions, providing promising directions for future research.