Anticoagulant peripheral blood and BM aspirate of SAA patients and healthy donors were obtained from Guangzhou First People’s Hospital. Patients were diagnosed with acquired AA according to the international AA Study Group Criteria (Marsh et al., 2009, Ball et al., 2009). Samples were centrifuged at 450 g for 5 min at 4°C; then plasma and BM supernatant were aliquoted and stored at −80°C until further processing.

Proteins and metabolites were extracted as previously described (Dunn et al., 2011, Broadhurst et al., 2011; Yan et al., 2020, Da et al., 2020). Briefly, 450 μl of methanol (Optima™ liquid chromatography–mass spectrometry (LC-MS) grade) was added to 150 μl of plasma or BM supernatant. The mixtures were thawed on ice with shaking for 3 min. The mixed solutions were centrifuged at 14,000 g for 10 min at 4°C. A quarter of each supernatant was subjected to centrifugal freeze-drying (Beckman Allegra 64R, Beckman Coulter, Brea, CA, USA) and then dissolved in 100 μl of MeOH:H₂O (v:v = 8:2; Pierce™ Water, LC-MS grade; Thermo Fisher Scientific, Waltham, MA, USA). After centrifugation at 14,000 g for 10 min at 4°C, the supernatant was collected for further detection.

The metabolites in plasma and BM samples were separated on a high-performance liquid chromatography (HPLC) system (Thermo Fisher Scientific, Waltham, MA, United States) with a C18 column (100 mm × 2.1 mm, 1.7 μm, Thermo Fisher Scientific) at a flow rate of 0.25 ml min⁻¹. Mobile phase A was 0.1% aqueous formic acid (v/v), and mobile phase B was methanol. Gradient elution was optimized as follows: 0–2 min, 2% B; 2–20 min, 2–100% B; 20–24 min, 100% B; 24–24.1 min, 100–2% B; and 24.1–30 min, 2% B. The column temperature was maintained at 35°C, and the injection volume was 5 µl.

The high-resolution MS for metabolite analysis was performed on a Thermo Scientific Q Exactive Plus mass spectrometry (Thermo Fisher Scientific, United States) with an electrospray ionization (ESI) source. Both negative and positive ion scan modes were carried out. The scan range was set at m/z 70–1,050 with an MS resolution of 70,000. The other MS parameters included the following: auxiliary gas flow, 10 arbitrary units; sheath gas flow, 35 arbitrary units; automatic gain control (AGC) target, 3e6 ions; capillary temperature, 320°C; S-lens radio frequency (RF), 50%; maximum injection time, 100 ms; and spray voltage, 3.5 kV.

Metabolites were identified by automatically comparing the precursor and product ions features (m/z) in the experimental samples with those of the chemical standards in the reference databases of Compound Discoverer 2.1.0.401 software (Thermo Fisher Scientific). The LC-MS/MS raw data were processed using Compound Discoverer 2.1 software to obtain the QC-normalized peak area. Then, the data were further processed in the MetaboAnalyst platform (https://www.metaboanalyst.ca/) including normalized by median and transformed by log(base 10) without scaling. Then, the normalized data were analyzed by partial least squares discriminant analysis (PLS-DA) and a volcano plot.

The samples for the proteomics study were processed and analyzed using the method and instrument parameters as described previously (Yan et al., 2020, Da et al., 2020), and proteins were identified using Proteomics Discovery Software (version 2.1, Thermo Fisher Scientific). Proteome data were analyzed in R (RStudio, Boston, MA, USA), and the pathway enrichment was analyzed using gseGO function in ClusterProfiler package (version 3.16.1) (Yu et al., 2012, Wang et al., 2012).

The original metabolome and proteome data can be found in the Supplementary Materials.

The publicly available single-cell RNA sequencing data (GSE145669) and microarray data (GSE3807) were obtained from the Gene Expression Omnibus (GEO). As for single-cell sequencing data quality control, cells that have more than 4,000 or less than 200 unique genes or have more than 25% mitochondrial genes were removed. Only genes expressed in 20 or more cells were used for further analyses. Then ribosome and mitochondria genes were filtered out. The pre-processed data matrix was transformed to a Seurat object using CreateSeuratObject function and normalized and log-transformed, and the log-transformed matrix was used for all downstream analyses in Seurat package (version 4.0.1). Principal component analysis (PCA) was used for dimensionality reduction using highly variable genes. With the use of the first fifteen principal components as input, the data in two dimensions were visualized with Uniform Manifold Approximation and Projection (UMAP). The gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA) were performed using the GSVA package (version 1.36.3) and ClusterProfiler package (version 3.16.1), respectively. Gene sets were downloaded from the MSigDB database or collected from document literatures. The differences in pathway enrichment scores between different clusters were calculated using LIMMA package (version 3.44.3). As for microarray data analysis, the expression matrix was downloaded by GEOquery (version 2.10), and gene annotation was performed to transform the probe name to symbol; then expression matrix was normalized with log2. The normalized matrix was analyzed for differential expression, and then heatmap was obtained by pheatmap package (version 1.0.12).

C57BL/6 and BALB/c mice were purchased from Hunan SJA Laboratory Animal Co., Ltd (China). Male C57BL/6 mice were crossed with female BALB/c mice to generate CByB6F1 mice. To induce immune-mediated BM failure, CByB6F1 mice aged 8–12 weeks received 5 × 10⁷ mixed lymph node cells and splenocytes from age- and sex-matched C57BL/6 mice or phosphate-buffered saline (PBS) intravenously. Mice were sacrificed at 3–4 weeks after cell transfer. All mice were housed and bred in Laboratory Animal Research Center, South China University of Technology (SCUT), under a specific pathogen-free condition.

BM mononuclear cells were isolated, and red blood cells were depleted using red blood lysis buffer (Beyotime Biotechnology, Shanghai, China; C3702). For surface markers staining, cells were incubated with purified anti-CD16/CD32 antibody (BioLegend, San Diego, CA, USA) for 15 min at 4°C and then stained with fluorochrome-conjugated antibodies for 20 min at 4°C. The antibodies against Ter119 (TER-119), CD45 (104), CD11b (M1/70), CD19 (6D5), NK1.1 (PK136), CD3 (17A2), CD8α (53–6.7), FasL (MFL3), and CD44 (IM7) were purchased from BioLegend; the antibodies against CD4 (GK1.5) and CD62L (MEL-14) were purchased from BD Biosciences (BD Biosciences, San Jose, CA, USA). For intracellular cytokine staining, cells were resuspended in RPMI-1640 with 10% fetal bovine serum and stimulated with Cell Stimulation Cocktail (plus protein transport inhibitors) (eBioscience, San Diego, CA, United States) at 37°C for 3 h. After surface marker staining, cells were fixed and permeabilized with Cytofix/Cytoperm™ Fixation/Permeabilization Kit (BD Biosciences) and then stained for anti-IFN-γ and anti-TNF-α (BioLegend). Flow cytometric analysis was performed on LSR Fortessa (BD Biosciences), and data were analyzed with the Flowjo software (BD Biosciences).

Data were analyzed using R (RStudio, Boston, MA, United States) and GraphPad Prism 8 (GraphPad Software, San Diego, CA, USA). Results in all figures were expressed as mean ± standard error of the mean. For the comparison of two groups, if they pass homogeneity of variance test (p > 0.05), statistics were analyzed using a two-tailed unpaired t-test with Welch’s correction; otherwise, a nonparametric test was used. p-Values less than 0.05 were considered significantly different: ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗p < 0.001.

Plasma proteomics has gained much interest for identifying possible prognostic markers in disease conditions. We collected the plasma from 15 healthy donors and 14 SAA patients and performed untargeted proteomics and metabolomics LC-MS analyses (Figure 1A) as previously reported (Want et al., 2006, O’Maille et al., 2006; Yuan et al., 2012, Breitkopf et al., 2012). The clinical characteristics of patients are summarized in Table 1. We mixed four to five plasma samples as one sample for proteomics analysis.

Plasma protein profile in SAA patients and HD is distinct and shows a good intra-group consistency by PCA and Pearson’s correlation analysis (Figure 1B; correlation coefficient >95%, Figure 1C). Of the 239 identified proteins, 49 had missingness and were removed from further analyses. Among the remaining proteins, 6 decreased and 15 increased proteins in SAA patients were observed (p < 0.05) (Figure 1D; Supplemental Material S1).

Upregulated proteins in SAA patients are enriched in immunoglobulin-mediated immune response (especially humoral immune response) and complement activation pathways by pathway enrichment analysis (Figure 1E). Moreover, complement molecules and serum amyloid A might participate in these pathways revealed by cnetplot analysis (Figure 1F).

Several significantly changed proteins between SAA patients and healthy donors, such as SHBG, APOA1, and PGK1, were closely associated with the metabolic state of patients (Figure 1D). Therefore, we profiled metabolomics of plasma from SAA patients and HD. We identified 225 known metabolites for further analyses. SAA patients and HD have distinct metabolite profiles, which were analyzed by PLS-DA (Figure 2A). A total of 50 metabolites were significantly changed (22.2% of 225 metabolites detected; Figures 2B,C); 33 metabolites were obviously increased, and 17 metabolites were significantly decreased (Figure 2C, Supplemental Material S2). Among them, drug-related metabolites (cefoperazone, norlidocaine, lidocaine N-oxide, lidocaine, metronidazole, omeprazole sulfone, cilastatin, resmethrin, and boldenone undecylenate) were significantly increased in plasma of SAA patients as compared with those of healthy donors showed by variable importance in projection (VIP) scores and a volcano plot (Figures 2B,C, Supplemental Material S2). Interestingly, sex hormones such as progesterone and testosterone isocaproate were also elevated in the plasma of SAA patients (Figure 2C, Supplemental Material S2). Notably, upregulated metabolites in SAA patients were enriched in steroidogenesis, arachidonic acid metabolism, and bile acid biosynthesis, while downregulated metabolites were enriched in the oxidation of branched-chain fatty acids, ketone body metabolism, beta oxidation of very-long-chain fatty acids, and caffeine metabolism pathways by metabolite set enrichment analysis (MSEA) (Figure 2D).

BM is the main target in SAA. Therefore, we also collected the BM supernatant of SAA patients (n = 5) to profile the metabolome. PLS-DA plot perfectly clustered BM supernatant and plasma (Figure 2E), showing that the metabolites in BM supernatant significantly differed from those in plasma. Drug-related metabolites (lidocaine, lidocaine N-oxide, and norlidocaine) were also increased in BM supernatant (Figures 2F,G, Supplemental Material S3). Interestingly, we found that 2-arachidonoylglycerol (2-AG) was substantially increased in both plasma and BM supernatant (Figure 2C,E,F) of SAA patients, suggesting that the endocannabinoid system may participate in the disease development. Upregulated metabolites in BM supernatant are enriched in bile acid biosynthesis, nicotinate, and nicotinamide metabolism and purine metabolism pathways revealed by MSEA. Meanwhile, metabolites associated with steroidogenesis and arachidonic acid metabolism pathways were downregulated, as well as galactose metabolism and caffeine metabolism in BM supernatant (Figure 2H).

To fully understand the functional role of these changes in the plasma of SAA patients, we combined the plasma proteomics and metabolomics data by joint pathway analysis (Basu et al., 2017, Duren et al., 2017) by analyzing upregulated proteins and metabolites. The results of the enriched 119 pathways were simultaneously plotted to screen pathways in terms of hypergeometric test p-values and pathway impact (Figure 3A). We found that the top ten pathways by p-value (top six) or impact (top four) were indicated the following pathways: glycolysis or gluconeogenesis, fructose, and mannose metabolism, purine metabolism, caffeine metabolism, nicotinate and nicotinamide metabolism, staphylococcus aureus infection, complement and coagulation cascades, HIF-1 signaling, cholesterol metabolism, and phenylalanine metabolism pathways (Figure 3A). Alteration in these pathways in SAA plasma provides clues that the coordination changes in proteins and metabolites may impact the system energy metabolism and immune response of SAA patients.

To find a possible link between systemic metabolic changes and CD8⁺ T-cell activation in AA patients, we further analyzed single-cell transcriptome data (GSE145669) of CD8⁺ T cells from AA patients and healthy controls (Zhu et al., 2021, Lian et al., 2021). Dimension reduction and unsupervised clustering classified CD8⁺ T cells into eight groups with distinct gene expression patterns (Figure 3B). Cells in clusters 0, 2, and 5 highly expressed the markers associated with the naïve T-cell phenotype, such as CCR7, SELL, and IL7R; cells in clusters 1, 3, 4, and 6 highly expressed the markers associated with memory or effector phenotypes, such as HLA-DRB1, KLRD1, NKG7, GZMA, and GZMB; and cells in cluster 7 highly expressed the genes associated with long-lived memory precursor cells (Kaech et al., 2003, Tan et al., 2003), such as KLRB1 and IL7R (Figure 3C).

Among the eight clusters, cells in clusters 1, 5, and 7 increased in AA patients (Figure 3D), which enriched genes involved in cell proliferation and survival, such as cell cycle phase, MYC target, and PI3K/AKT/MTOR pathways. They also express genes associated with proinflammation pathways that are involved in AA, such as TNFα signaling (Figure 3E). These results suggest that CD8⁺ T cells in AA patients were different from healthy controls in both phenotype and intrinsic pathways, which may be responsible for the systemic metabolic changes in AA.

We further discovered the elevated expression of T-cell activation (CD38, GZMB, GZMA, CD27, and NKG7) genes in BM T cells from AA patients by analyzing microarray data (GSE3807) (Franzke et al., 2006, Geffers et al., 2006) (Figure 4A). Previous reports showed that CD38 drives T-cell proliferation, cytokine secretion, and activation (Muñoz et al., 2008, Mittelbrunn et al., 2008). Also, CD38 is closely related to the metabolic state and regulates T cells' fate due to its multifunctional enzymatic activity such as NAD/NMNase (Hogan et al., 2019, Chini et al., 2019).

On this basis, we hypothesized that CD38 might be involved in the pathogenic process of CD8⁺ T cells in SAA. We defined CD38⁺ cells by the expression value of CD38 > 0 in the single-cell transcriptome data (GSE145669) and observed that CD38⁺CD8⁺ T cells are mainly grouped within clusters 1, 5, and 7 (Figure 4B). Moreover, we found that CD38⁺CD8⁺ T cells were significantly increased in patients (Figure 4C).

Further, we compared the pathway variation between CD38⁺CD8⁺ T cells and CD38⁻CD8⁺ T cells within clusters 1, 5, and 7. CD38⁺CD8⁺ T cells express genes enriched with the pathway participated in SAA by GSVA, such as TNF-α signaling, indicating that CD38⁺CD8⁺ T cells might be important in SAA (Figure 4D). CD38⁺CD8⁺ T cells also express hallmark genes involved in hypoxia, glycolysis, fatty acid metabolism, and complement pathways (Figure 4D). Moreover, we observed that T-cell microarray data displayed increased hypoxia-associated gene expression (CD46, CITED2, and HIGD1A) in AA patients (Figure 4A), which further demonstrated that CD38⁺CD8⁺ T cells are activated and undergo metabolic changes in AA. These data indicated that CD38⁺CD8⁺ T cells might be important in disease development.

To demonstrate the pathogenic role of CD38⁺CD8⁺ T cells in AA, we then constructed a murine immune-mediated AA model by infusing C57BL/6 lymphocytes into CByB6F1 mice, which induces AA symptoms such as hypocellular BM and lymphocyte infiltration in the BM (Krenger et al., 2000, Rossi et al., 2000). We observed a significant decline in the total number of BM cells (Figure 5A), as well as a decreased percentage of Ter119⁺ cells (Figure 5B) and B cells (Figure 5C). However, the percentage of T cells (Figure 5D) was significantly increased, and CD8⁺ T cells were predominant in the BM (Figure 5E).

Consistent with single-cell transcriptome data, CD38⁺CD8⁺ T cells significantly increased in AA mice, as well as the expression and MFI of CD38 (Figure 5F). Moreover, most of CD38⁺CD8⁺ T cells were CD62L⁻CD44⁺ (Figure 5G), especially in AA mice. It has been evidenced that CD8⁺ T cells usually achieve hematopoietic cell injury through the Fas/FasL pathway (Young et al., 2006, Calado et al., 2006). Hence, we analyzed the level of FasL on CD38⁺CD8⁺ T cells between AA and CTRL. As expected, the expression of FasL was significantly higher in AA mice, which was consistent with a previous report (Figure 5H). Also, BM CD8⁺ T cells from AA mice have significantly elevated IFN-γ secretion ability (Figure 5I). We also detected the secretion of IFN-γ and TNF-α by CD38⁺CD8⁺ T cells and CD38⁻CD8⁺ T cells in the BM of the induced AA murine model. We found that CD38⁺CD8⁺-derived T cells secreted more IFN-γ and TNF-α, indicating that CD38⁺CD8⁺ T cells might be the more potent pathogenic immune cells during AA development (Figure 5J). These data suggest an activated phenotype of CD38⁺CD8⁺ T cells in the BM of AA.