Three datasets of HSC expression profiles from young and aged mice (GSE100906, GSE70657 and GSE100426, dataset information is shown in Supplementary Table S1) were downloaded from the GEO database. To minimize batch effect between datasets, we integrated three datasets following the instructions of Seurat (Butler et al., 2018; Stuart et al., 2019). Briefly, the FindVariableFeatures function was used to identify the top highly variable 2000 genes and the IntegrateData function was used to move the batch effect and correct the datasets (dims = 1:30). Then the ScaleData function of Seurat was used to calculate the scaling expression values and the RunPCA function (npcs = 30) was used to perform a PCA dimensionality reduction by Uniform Manifold Approximation and Projection (UMAP). The FindNeighbors and FindClusters Seurat functions were used to cluster the cells (resolution = 0.8) (See Supplementary Method for details).

DEG identification was performed using the package DESeq2 (Love et al., 2014). |Foldchange| > 1.5 and p < 0.05 were set as the cutoffs to screen for DEGs. Firstly, bone marrow cells were isolated from young (4∼6 weeks, n = 3, C57BL/6) and aged (18∼20 months, n = 3, C57BL/6) mouse femora and tibiae by flushing bones with 1–3 ml phosphate-buffered saline (PBS) (without Mg²⁺ and Ca²⁺) supplemented with 2% fetal bovine serum. Then the bone marrow cells were resuspended in cold ammonium chloride solution and incubated on ice for 10 min to lyse red blood cells (RBCs). Following RBC lysis, the cells were washed in PBS and then incubated with flow cytometry antibodies (Lineage-AF700, c-kit-PECY5, Sca-1-FITC, CD150-PE, CD48-PECY7) for 30 min. LT-HSCs (Lineage- Sca1+ c-Kit + CD150 + CD48-) were collected with an SH800S Cell Sorter. As the number of sorted LT-HSCs was small, we amplified the transcriptome in single cells by Single Cell Sequence Specific Amplification Kit (P621-01/02, Vazyme, Nanjing, China) and then detected the gene expression by a fluorescent quantitative PCR according to the manufacturer’s instructions. Firstly, the amplification primers of DEGs were mixed together to prepare an Assay Pool (primer concentration 0.1 µM). Then single LT-HSC were transferred into PCR strips with Assay Pool, RT/Taq enzyme and Reaction Mix. After 5 min of freezing at −80°C, the PCR strips were kept at 50°C for about 60 min to perform reverse transcription. Then the cDNA was denatured at 95°C for 15 s, annealed at 60°C for 15 min, and 15 sequence-specific amplification cycles were performed. After pre-amplification of transcriptome, qRT-PCR was performed following the one-step protocol of the SYBR Green master mix (Vazyme) in a Real-time system according to the manufacturer’s instructions (95 C*15 s, 60 C*30 s, 40 cycles). The primer sets for the detection of DEGs are listed in Supplementary Table S3.

Functional enrichment analysis of DEGs was performed with DAVID website (Huang et al., 2008; Huang et al., 2009). Gene Ontology (GO) term categories included biological process (BP), molecular function (MF) and cellular component (CC). STRING website was used to predict the interactions among DEGs (Szklarczyk et al., 2019). Cytoscape software was applied to establish the PPI network (Shannon et al., 2003) and MCODE was applied to identify the hub genes (Bader and Hogue, 2003). CytoHubba, a plug in Cytoscape, was applied to explore important nodes/hubs in the interactome network (Chin et al., 2014).

The cell cycle phase (G1, S, or G2M) of each cell was identified by calculating the cell cycle phase score based on canonical markers using the Seurat package (version 3.1) (Butler et al., 2018; Stuart et al., 2019). The canonical cell cycle gene set was defined with “cell cycle process” from the GO annotation and identified as cycling in HeLa cells (Whitfield et al., 2002) and it can be downloaded from GSEA MSigDB version 7.1 (https://www.gsea-msigdb.org/gsea/msigdb) (Subramanian et al., 2005; Liberzon et al., 2015).

The single cell transcriptome profiles of young and aged mouse bone marrow based on the microfluidic droplet platform (10x Genomics) were shared by the Tabula Muris Consortium (tabula-muris-senis.ds.czbiohub.org) (Schaum et al., 2018; Almanzar et al., 2020) and can be downloaded from the Gene Expression Omnibus database (GSE109774 and GSE132042). CellphoneDB software was used for cellular interaction analysis (Vento-Tormo et al., 2018; Efremova et al., 2020). Significant ligand–receptor pairs were filtered with p < 0.05. To identify which cell population accounted for the difference in inflammatory cytokine levels, CellChat software was applied to analyze outgoing communication patterns of secreting cells (Jin et al., 2021). IdentifyOverExpressedGenes (thresh.p = 0.05) and identifyOverExpressedInteractions functions were used to identify over-expressed ligands or receptors and over-expressed ligand-receptor interactions. ComputeCommunProb function (methods for computing the average gene expression per cell group was set as “triMean”) was used to infer the biologically significant cell-cell communication (See Supplementary Method for details).

LEGENDPlex assays (BioLegend) were performed to detect the inflammatory cytokine and chemokine levels of bone marrow from young (4∼6 weeks, n = 3, C57BL/6) and aged (18∼20 months, n = 3, C57BL/6) mice. The femurs and tibias were cut at both ends, and the bone marrow was flushed out with 1 ml PBS. After 5 min of centrifugation at 400 g, the supernatant was collected. Samples were detected according to the manufacturer’s instructions by flow cytometry.

Colony formation unit assays were performed as follows. One hundred FACS-isolated HSCs (Lineage-Sca1+c-Kit+) from young (4∼6 weeks, n = 3, C57BL/6) and aged (18∼20 months, n = 3, C57BL/6) mice were plated in methylcellulose (M3434 for BFU-E, CFU-GM, CFU-GEMM measurement and M3630 for CFU-B measurement, Stem Cell Technologies) and enumerated on day 10. To test the effect of inflammatory cytokines on hematopoietic cell lineage, either DMSO vehicle or IL-1β (25 ng/ml), TNF-α (1 μg/ml), AS101 (2 μg/ml, IL-1β inhibitor) and R7050 (2 μg/ml, TNF-α receptor antagonist) were added to the methylcellulose media.

To characterize the differences in HSC subpopulations, three datasets (Supplementary Table S1) of HSC single cell expression profiles from young and aged mice were integrated together using Seurat (Butler et al., 2018; Stuart et al., 2019) to correct batch effect. A total of 5 clusters were identified (Figures 1A,B) by using UMAP based on their markers and gene set enrichment analysis (Figures 1C,E). To assess the impact of aging on HSC population, we compared the proportion of different clusters in young and aged HSCs (Figure 1D). Interestingly, Cluster 4, a proliferation-associated cluster marked by Cdca7, decreased significantly upon aging. Besides, Cluster 2 (an apoptosis-associated cluster marked by Stat3) and Cluster 5 (an interferon-associated cluster marked by Ifngr1) were increased upon aging. Ifngr1 encodes the ligand-binding chain of the gamma interferon receptor and IFN-γ has been reported to selectively promote differentiation of myeloid-biased HSCs (Matatall et al., 2014). As a whole, these results highlighted an exhaustion of HSCs being able to give rise to HSC itself without differentiation, which is a hallmark of HSC aging and indicated that Stat3 and Ifngr1 were two markers for accumulation of inflammatory and apoptosis-biased state HSCs.

In order to elucidate the overall alterations of HSC aging, we next analyzed HSCs as a whole. According to the cut-off criterion (|Foldchange| > 1.5 and p < 0.05), 453, 614 and 297 DEGs between young and aged HSCs were identified from GSE100906, GSE70657 and GSE100426, respectively (Figure 2A). Common DEGs were defined as genes that were significantly upregulated or downregulated in at least two datasets. By employing integrated bioinformatics analysis, 56 common upregulated genes and 51 common downregulated genes were identified (Figure 2B, Supplementary Table S2). For instance, Birc5 and Kpna2 were downregulated, while Clu, Selp and Sdpr were upregulated in HSCs during aging. To confirm the results of common DEGs, the relative expression levels of top 30 DEGs were analyzed by qRT-PCR (Figure 2C). We found that the PCR results for approximately 75% of the genes were consistent with our bioinformatics analysis (p < 0.05).

To further delineate the functional changes that occur during HSC aging, functional enrichment of the common DEGs was performed by using the DAVID gene annotation tool (Figure 3A). For BP, the upregulated genes were mainly enriched in transcription, DNA-templated and regulation of transcription from RNA polymerase II promoter, and the downregulated genes were predominantly enriched in cell cycle, mitotic nuclear division, cell division and protein phosphorylation. This was consistent with the previous findings that cell cycle-related genes dominated the transcriptomic variability of aging and that aged HSCs underwent fewer cell divisions than young HSCs (Janzen et al., 2006; Nygren and Bryder, 2008; Kowalczyk et al., 2015). KEGG pathway enrichment analyses revealed that the upregulated DEGs were mainly involved in osteoclast differentiation and TNF signaling pathway, and the downregulated DEGs were mainly involved in cell cycle, oocyte meiosis and p53 signaling pathway. p53 is implicated in regulating HSC aging and quiescence (Dumble et al., 2007) and regulating p53 can help to maintain hematopoietic cells during oxidative stress (Jung et al., 2013). These results highlighted the loss of self-renewal ability and quiescence in aged HSCs.

To better understand the interactions among the common DEGs, a PPI network with 67 nodes (27 upregulated and 40 downregulated) and 413 edges was generated with the STRING tool (Figure 3B). Two significant modules were also identified based on the degree of importance in Cytotype MCODE. Module 1 contained 25 nodes and 286 edges (Figure 3C), and the expression of all the nodes was downregulated in aged HSCs. KEGG pathway enrichment analyses revealed that the DEGs in Module 1 were mainly enriched in the cell cycle. Module 2 contained 11 nodes and 17 edges (Figure 3D), and the DEGs in Module 2 were mainly enriched in cell adhesion molecules and the hematopoietic cell lineage. In addition, highly connected nodes in the network are likely to have significant functional importance and were defined as hub genes. By utilizing CytoHubba in Cytotype, we identified the top five upregulated hub genes (Jun, Aldh1a1, Egr1, Cd38 and Junb) and the top five downregulated hub genes (Aurka, Ccna2, Ccnb2, Cdk1 and Birc5). The alterations and potential functions of the hub genes were summarized in Table 1.

As the enrichment analysis of common DEGs and PPI network analysis suggested that the cell cycle was strongly associated with functional decline in aged HSCs, we next compared the expression levels of cell cycle-associated genes. Among 107 DEGs, 24 genes were associated with the cell cycle and most of them (22 genes), including Chek2, Kif20b and Clasp2, were downregulated in at least 2 datasets (Figure 4A). Subsequently, the cell cycle phase of young and aged HSCs was identified by calculating the cell cycle phase score based on canonical markers in both young and aged mice (Figure 4B and Supplementary Figure S1A). The frequency of cells in the G1 cluster significantly decreased in aged versus young cells (13.5 vs. 6.8%; p < 0.05, shown in Figure 4C).

To further identify the signaling pathways that drove HSCs to transit through G1 phase quickly, we evaluated genetic alterations in a cell cycle pathway diagram that was publicly available at https://www.genome.jp/kegg (Kanehisa and Goto, 2000). We found that the Skp2-induced signaling pathway (Skp2→Cip1→CycA/CDK2→DP-1) was significantly downregulated during the aging process; this change promoted the synthesis of mRNAs and proteins that are required for DNA synthesis during S phase (Figure 4D).

The PPI network analysis confirmed significant changes of adhesion molecules, highlighting an important role of bone marrow microenvironment in the aging process. Therefore, we reanalyzed two single cell transcriptome profiles of young and aged mouse bone marrow (GSE109774 and GSE132042) and investigated how the cellular composition of the bone marrow changed with age. The percentages of erythroblasts, granulocytopoietic cells and granulocytes increased with age, while the percentages of macrophages, late pro-B cells, monocytes and T cells decreased with age (Figure 5A). Cell-cell communication mediated by receptor-ligand complexes is important for stem cell biological processes (Efremova et al., 2020). To assess alterations in intercellular communication during aging, cellphoneDB software was used to predict interactions among different cell types from the single-cell RNA-seq data (Vento-Tormo et al., 2018; Efremova et al., 2020). Significant predicted interactions were assessed separately for young and aged mice (Figure 5B), and the differences were used to infer alterations in intercellular interactions (Figure 5C). The predicted interactions with hematopoietic stem/progenitor cells (HSPCs) were highest in monocytes and lowest in erythroblasts in both young and aged mice. The comparison between the young and aged groups indicated a significant increase in intercellular communication between HSPCs and proerythroblasts, and a decrease in intercellular communication between HSPCs and monocytes, granulocytopoietic cells and granulocytes (Figure 5C).

As the surrounding cells of HSPCs are potential niche components, the ligands expressed by these populations and the receptors expressed by HSPCs were considered in the downstream analysis. Among a total of 68 heterologous ligand–receptor pairs detected between HSPCs and other cell types, immune response, inflammation response and signal transduction were the predominant biological processes involved (Figure 5D). In addition, ligand–receptor pairs were involved in some canonical signaling pathways (such as the TNF signaling pathway, NF-kappa B signaling pathway and PI3K-Akt signaling pathway) and certain aging-associated diseases (such as Alzheimer’s disease).

The ligand–receptor gene pairs exhibited differential expression patterns between young and aged populations when coupled with HSPCs. Interactions with HSPCs changed significantly in proerythroblasts, monocytes, granulocytopoietic cells and granulocytes, and the significant differentially expressed ligand–receptor gene pairs of these four groups of cells are shown in Figure 5E. Adhesion complexes, such as aMb2 complex-Icam1 and aLb2 complex-Icam1, were differentially expressed during HSC aging. Moreover, inflammatory ligand-receptor pairs, such as CXCL2-CXCR2, IL15-IL15 receptor, CCL2-CCR2 and CCL5-CCR5, were upregulated during HSC aging.

The upregulated inflammatory ligand-receptor pairs in aged HSCs suggested increasing inflammation levels in the bone marrow niche during aging. To test our hypothesis, a Cytoplex Assay was performed to measure the levels of inflammatory cytokines and chemokines in the bone marrow. Most inflammatory cytokines and chemokines, including TNF-α, IFN-β, IFN-γ, IL-1α, IL-1β, IL-6, IL-17α, IL-23, CCL2, CCL4 and CXCL1, were upregulated in aged bone marrow niche, while CCL20 and CXCL10 were downregulated in the aged bone marrow niche (Figure 6A). To assess the effect of inflammatory cytokines on cellular senescence, HSCs were cultured in methylcellulose with either inflammatory cytokines or their inhibitors. TNF-α and IL-1β promoted myeloid-biased differentiation and inflammatory pathway blockade may rejuvenate aged HSC functions and increase B cell output (Figure 6B).

To further determine which cell population accounted for the difference in inflammatory cytokine levels, the outgoing communication patterns of secreting cells were analyzed (Figure 6C). The secretion of HSPCs, erythroblasts and basophils contributed to high CCL levels in the aged bone marrow niche, while the secretion of T cells, promonocytes, macrophages, erythroblasts, and basophils contributed to high IFN-γ levels. In addition, the secretion of some interleukins and complement factors was upregulated in proerythroblasts and immature B cells. Interestingly, the levels of inflammatory cytokines secreted by HSPCs, including CCL and IL2, were increased, highlighting the important role of autocrine signaling during aging process. Circle plots of the cell-cell communication network also revealed that inflammatory pathway networks (such as IL2, CCL and complement; Figure 6D, Supplementary Figure S1B) became more interconnected during aging, with some ligand-receptor communication enhanced significantly (such as IL2-IL2RB/IL2RG and CCL5-CCR1; Supplementary Figure S1C).