BM aspirates and PB samples were collected from 7 adult AA patients at diagnosis (pre-ATG; AA^PRE) and 7 age- and sex-matched healthy donors (HDs). In addition, BM aspirates from 3 of 7 AA patients were collected six months after start of first-line IST (post-ATG; AA^POST), which consisted of 4 days 40mg/kg ATGAM^® (hATG; Pfizer) and 5mg/kg/day ciclosporin for at least six months. BM and PB collection was performed after given informed consent in accordance with the Declaration of Helsinki. Clinical characteristics of patients and HDs are outlined in Table 1 . All patients were treatment-naïve at diagnosis, had an indication for treatment, did not have an additional underlying disease and were diagnosed and treated according to Dutch guidelines (3). Patients under the age of 40 were screened for Fanconi anemia to rule out a congenital form of BM failure. 1 patient was classified as non-severe AA (NSAA), 2 patients as severe AA (SAA) and 4 patients as very severe AA (VSAA), defined as described in international guidelines (2). Complete response to IST was defined as complete normalization of PB counts. Partial response was defined as transfusion independence and a neutrophil count of at least 0.5·10⁹ neutrophils/L. This study was approved by the Medical Ethical Committee of the Leiden University Medical Center (protocol number B20.037).

Bone marrow mononuclear cells (BMMCs) or PB mononuclear cells (PBMCs) were isolated from BM aspirates or PB, stained with a panel of 39 validated (12–14) metal-isotope tagged antibodies and measured as described in Supplemental Methods . For each acquired mass cytometry data file, single, live CD45⁺ cells were gated in FlowJo version 10.7.1 (BD Biosciences) as shown in Supplementary Figure 1 . A median of 0.5·10⁶ (range 0.4-0.9·10⁶) cells was gated for each AA^PRE sample. For every HD or AA^POST sample, a median of 0.5·10⁶ (range 0.3-0.9·10⁶) or 0.3·10⁶ (range 0.3-0.6·10⁶) cells was gated respectively.

The single, live CD45⁺ cells were studied using two analyses. In the first analysis, BMMC samples from all AA^PRE and HDs were studied. 7.7·10⁶ single, live CD45⁺ cells were sample-tagged, hyperbolic ArcSinh transformed with a cofactor of 5 and imported in Cytosplore (15) for dimensionality reduction. A five-level hierarchical stochastic neighbor embedding (HSNE) analysis was performed with default perplexity and iterations (30 and 1000, respectively; see Supplementary Figure 2 for outline of HSNE analysis pipeline). Three of 39 markers in the mass cytometry panel (HLA-DR, CD5 and KLRG1) were not used in the HSNE clustering analysis due to small variations in marker distribution between experiments (data not shown). Major immune lineages were identified at the overview level of the HSNE analysis and were defined as: CD4⁺CD3⁺CD7⁺ CD4⁺ T-cells, CD8a/b⁺CD3⁺CD7⁺ CD8⁺ T-cells, CD7^-CD20/IgM⁺ B-cells, CD3⁺CD7⁺TCRgd⁺ non-conventional T-cells (NCTs), CD3^-CD7⁺ innate lymphoid cells (ILCs) including NK-cells and Lin^-cKit⁺ HSPCs. All remaining cells were considered myeloid cells. Subsequently, each major immune lineage was selected for detailed analysis at the data level of the HSNE analysis using t-distributed stochastic neighbor embedding (tSNE) analyses with up to 0.5·10⁶ landmarks. Clusters with distinct surface marker expression were identified within each major immune lineage by Gaussian mean shift (GMS) clustering at the data level. Each cluster consisted of at least 1000 cells ( Supplementary Figure 2 ).

In the second analysis, paired BMMC samples from 3 AA^PRE and AA^POST were analyzed using a comparable approach, for detailed description see Supplemental Methods .

After clustering in both analyses, .fcs files of major immune lineages and clusters were exported from Cytosplore and loaded into R version 4.0.5 for data visualization and analysis. Cytofast (16) was used to generate heatmaps and to compare lineage or cluster frequencies between groups. Principal component analyses (PCA) and a Spearman’s rank correlation analysis were performed using the packages factoextra version 1.0.7 (17) and corrplot version 0.92 (18), respectively.

Mann-Whitney U tests performed in R version 4.0.5 were used to compare immune cell frequencies between two groups. P-values ≤0.05 were considered statistically significant. Since this was an exploratory study, no adjustments for multiple testing were applied. *P ≤ 0.05; **P ≤ 0.01; ***P ≤ 0.001.

First, we studied the immune cell composition in the BM of the 7 AA^PRE and 7 age- and sex-matched HDs. A five-level HSNE analysis was performed to visualize all 7.7·10⁶ CD45⁺ cells (median 0.5·10⁶ CD45⁺ cells per sample). At the overview level of the HSNE analysis, all major immune lineages and HSPCs could be identified by clustering based on lineage markers and cKit ( Figure 1A ). All AA^PRE and HD samples were homogeneously distributed per group over the HSNE analysis ( Supplementary Figure 3 ). Visualization of the AA^PRE and HDs in the overview level already revealed clear differences in the immune cell composition.

An unsupervised PCA confirmed a clear separation in the immune cell composition of AA^PRE and HDs ( Figure 1B ). Next, the immune cell composition was studied in further detail at the major immune lineage level. Since we aimed to uncover how the immune cell composition differs in BM from AA^PRE and HDs and total cell counts are decreased in the BM in AA, we studied relative cell frequencies instead of absolute cell counts. Major immune lineage frequencies of CD45⁺ cells were determined ( Figure 1C ). HSPCs could not be distinguished enough at this level of the HSNE analysis and were therefore merged with myeloid cells. In agreement with disease characteristics, myeloid cells were scarce in AA^PRE compared to HDs (median 0.6% versus 33% of CD45⁺ cells). B-cells were also reduced in AA^PRE compared to HDs (median 9% versus 16% of CD45⁺ cells). In contrast, CD4⁺ and CD8⁺ T-cells were increased in AA^PRE compared to HDs, which could be due to the relative decrease in myeloid and B-cells. Finally, non-conventional T-cell (NCT) and innate lymphoid cell (ILC) frequencies did not significantly differ between AA^PRE and HDs.

Since myeloid, B-cell, CD4⁺ T-cell and CD8⁺ T-cell frequencies were significantly different in BM from AA^PRE compared to HDs, we studied each of these lineages at the immune subpopulation level. First, the myeloid cells including HSPCs were selected and analyzed at the four successive levels of the HSNE analysis ( Figure 2 , Supplementary Figure 2 ). In this manner, multiple HSPC and myeloid cell populations could be distinguished based on the expression of defining cell surface markers ( Figure 2A , right panels). Twenty phenotypically distinct clusters were identified ( Figure 2B ), 8 of which were significantly different in AA^PRE compared to HDs, 1 had a p-value of 0.053 and 11 clusters did not significantly differ between both groups. Analysis of the 8 significant clusters and cluster 9 showed a strong reduction in all cKit^+/dim clusters (#9-11) in AA^PRE, confirming HSPCs depletion in AA^PRE ( Figure 2C ). Furthermore, we observed a decrease of myeloid cells compatible with CD123⁺ plasmacytoid dendritic cells (#13), CD163⁺ macrophages (#16), CD11c⁺ conventional dendritic cells (#17) and CD14⁺CD16^- classical monocytes (#18) in AA^PRE. The most striking observation was a substantial increase in CD16⁺ myeloid cells in AA^PRE, which lacked expression of the neutrophil marker CD15 and represented a median of 41% (range 16-83%) of all myeloid cells ( Supplementary Figures 4A, B #14;19-20 combined; arrows in Figure 2A ). Two CD16⁺ myeloid cell clusters (#19-20) resembling CD14^-CD16⁺ non-classical and CD14⁺CD16⁺ intermediate monocytes were significantly increased in AA^PRE ( Figure 2D ). In contrast, CD16⁺ myeloid cells only comprised a median of 5% (range 2-19%) of all myeloid cells in HDs ( Supplementary Figures 4A, B ).

Detailed analyses of all B-cells revealed phenotypical heterogeneity within the B-cell compartment ( Figure 3A ). 20 B-cell clusters were identified ( Supplementary Figure 4C ). CD20⁺IgM^-CD127⁺ B-cells (#2) and CD20^-IgM^- B-cells expressing variable levels of CD127 (#3-4), compatible with a group of early B-cells (19, 20), were markedly reduced in AA^PRE compared to HDs ( Figure 3B ). We also observed a decrease in B-cells with absent or low expression of the chemokine receptor CCR6 ( Figure 3B , #14;17-20). Strikingly, B-cells expressing high levels of CCR6 were strongly increased in AA^PRE compared to HDs ( Figure 3C , #12-13;16). Combined, CCR6⁺⁺ B-cells constituted a median of 64% (range 38-86%) of all B-cells in AA^PRE and only a median of 10% (range 1-26%) of all B-cells in HDs ( Supplementary Figures 4C, D , #5;9;12-13;16 combined). Clusters 12 and 16 co-expressed CD27 in line with previous descriptions of memory B-cells (21–23) ( Supplementary Figure 4C ). In contrast, cluster 13 lacked expression of CD27, compatible with a population of immature or mature B-cells (24). In addition, the CCR6⁺⁺ B-cells expressed different levels of IgM, indicating that these cells were in multiple stages of development.

Hierarchical clustering of all CD4⁺ T-cells classified 35 phenotypically distinct CD4⁺ T-cell clusters in different stages of development ( Supplementary Figures 4E–G ). No significant differences were observed in naïve (CD45RA⁺CCR7⁺) CD4⁺ T-cell clusters. In contrast, multiple central memory (CM; CD45RO⁺CCR7⁺) and effector memory (EM; CD45RO⁺CCR7^-) CD4⁺ T-cell clusters were significantly reduced in AA^PRE ( Figure 4A ). The most remarkable difference was an increase in CM and EM CD4⁺ T-cells with a Th17-like phenotype based on the co-expression of CCR6 and/or CD161 (25), suggestive of Th17-skewing in AA^PRE. Six CCR6⁺CD4⁺ T-cell clusters were significantly increased in AA^PRE compared to HDs ( Figure 4B ). Together, the CCR6⁺CD4⁺ T-cells represented a median of 11% (range 5-18%) of CD4⁺ T-cells in AA^PRE, compared to a median of 2% (range 2-6%) of CD4⁺ T-cells in HDs ( Supplementary Figures 4F, G , #13;16;18;20;22;34-35 combined).

Hierarchical clustering of all CD8⁺ T-cells identified 45 phenotypically distinct CD8⁺ T-cell clusters ( Supplementary Figures 4H–J ). CM and EM CD8⁺ T-cells were predominantly reduced in AA^PRE ( Figure 4C , Supplementary Figures 4I, J ). A cluster of CD8⁺ T-cells (#36) resembling a population of naïve but activated (CD38⁺) T-cells was significantly increased in AA^PRE ( Figure 4D ). The most striking observation within the CD8⁺ T-cell compartment in AA^PRE was an increase of EMRA CD8⁺ T-cells positive for KLRG1 resembling highly differentiated cytotoxic T-cells (26). These KLRG1⁺EMRA CD8⁺ T-cells constituted a median of 29% (range 15-45%) of all CD8⁺ T-cells in AA^PRE, but only a median of 9% (range 7-33%) of all CD8⁺ T-cells in the HDs (#10-22;24-26;28-29 combined; Supplementary Figures 4I–J ). Five KLRG1⁺EMRA CD8⁺ T-cell clusters (#12;16-17;19;28) were significantly elevated in AA^PRE compared to HDs ( Figure 4D ).

To capture the relationships across the immune cell clusters that characterize AA^PRE, we performed a Spearman’s rank correlation analysis on the myeloid, B-cell and T-cell clusters with significant differential presence between AA^PRE and HDs. A matrix of Spearman’s correlation coefficients visualized two large networks of immune cells ( Figure 5A ). In line with the analysis of the major immune lineages above, network 1 showed a correlation between CD16⁺ myeloid cells, CCR6⁺⁺ B-cells, CCR6⁺ memory CD4⁺ T-cells, CD45RA⁺CCR7⁺CD38⁺ CD8⁺ T-cells and KLRG1⁺EMRA CD8⁺ T-cells, and was significantly more abundant in AA^PRE compared to HDs (median 22% versus 5% of CD45⁺ cells, Figure 5B ). In contrast, network 2 was significantly reduced in AA^PRE compared to HDs (median 6% versus 23% of CD45⁺ cells). Analysis of paired PB and BM samples from all 7 AA^PRE demonstrated that all subpopulations of network 1 are detectable in PB ( Supplementary Figure 5 ).

Since IST can result in recovery of hematopoiesis, we investigated the effect of IST on the immune cell composition in BM of 3 AA patients pre- and six months post-ATG. All patients still received ciclosporin and were transfusion-independent six months post-ATG. All patients had improvement of PB counts and eventually responded to IST, however, only one patient (AA3) showed complete response six months post-ATG (8.3 mmol hemoglobulin/L, 149·10⁹ platelets/L and 1.8·10⁹ neutrophils/L). One patient (AA2) had major partial response six months post-ATG (5.8 mmol hemoglobulin/L, 87·10⁹ platelets/L and 1.9·10⁹ neutrophils/L), while the third patient (AA1) showed minor partial response six months post-ATG (5.4 mmol hemoglobulin/L, 18·10⁹ platelets/L and 0.3·10⁹ neutrophils/L).

A total of 3.0·10⁶ CD45⁺ cells (median 0.5·10⁶ or 0.3·10⁶ cells per AA^PRE or AA^POST, respectively) were analyzed using a four-level HSNE analysis. At the overview level, the global immune cell composition in AA^PRE and AA^POST was visualized and all major immune lineages were identified ( Figure 6A ). HSPCs were not observed, suggesting that the frequency of HSPCs was too low to be detected. Unsupervised PCA presented a clear separation of AA^PRE and AA^POST ( Figure 6B ), indicative of a considerable shift in the immune cell composition in the BM six months post-ATG. Analyses at the major immune lineage level demonstrated that myeloid cells were strongly increased or normalized in all AA^POST ( Figure 6C ). In addition, CD4⁺ T-cell were normalized, while B-cells were still below normal levels in all AA^POST. CD8⁺ T-cells remained high in the patient (AA1) with the lowest blood counts six months post-ATG ( Figure 6C ).

To analyze the effect of IST on the disease-specific immune cell network identified in AA^PRE (network 1), we studied the presence of its corresponding immune cell subtypes in AA^POST. Cluster frequencies were compared to frequencies observed in the HDs. Strikingly, this revealed that CD16⁺ myeloid cell frequencies normalized in all AA^POST and CD45RA⁺CCR7⁺CD38⁺ CD8⁺ T-cell frequencies decreased in all AA^POST ( Figure 6D ). CCR6⁺⁺ B-cell, CCR6⁺ memory CD4⁺ T-cell and KLRG1⁺EMRA CD8⁺ T-cell frequencies showed a trend toward normalization in the two patients (AA2 and AA3) with complete response or major partial response six months post-ATG. In contrast, these clusters were still frequent or increased in the patient (AA1) with minor partial response six months post-ATG ( Figure 6D ).