A single-cell time series dataset capturing the dynamic inflammatory response of HPSCs

Biological responses such as acute inflammation are dynamic processes in which cells, tissues, and organisms undergo different phases of sensing differences, responding to these changes, and recovering upon successful response. Yet, often only single time points of these responses are investigated. Quiescent HSCs respond to inflammation by increased proliferation, which can be mimicked by treating mice with single pro-inflammatory cytokines, such as IFNα. To gain a better understanding of the dynamics of the HSCs response to acute inflammation, we performed a time series experiment to cover the sensing, response, and recovery phases of the inflammatory response of HSCs to treatment with IFNα. Whereas at 3 h post-treatment, HSCs showed the first signs of sensing IFNα by increasing expression of interferon-stimulated genes (ISGs) (Fig 1A), only at 24 h posttreatment, HSCs reached a peak in increased proliferation (Fig 1B) (Essers et al, 2009; Pietras et al, 2014). 72 h posttreatment, ISG expression was back to baseline (Fig 1A) and most of the HSCs returned to quiescence (Fig 1B). Unfortunately, inflammation does not only lead to increased proliferation of HSCs but is also accompanied by increased expression of several cell surface protein markers used to identify different cell types within the HSPC compartment. The most well-known example is the increase in the stem cell marker Sca-1 (Essers et al, 2009; Pietras et al, 2014; Kanayama et al, 2020). Using conventional marker-based flow cytometry including Sca-1, an increase in the more stem-like LSK (Lin⁻ Sca-1⁺ cKit⁺) cells and a decrease in the more committed LS⁻K (Lin⁻ Sca-1⁻ cKit⁺) progenitors is observed in response to inflammation (Fig 1C–E). However, because of the increase in protein expression of Sca-1 (Fig 1F and G) it is hard to predict whether these changes in frequency reflect an actual increase/decrease in cell frequency or are the result of contamination because of changes in protein expression of the cell markers in a given population or even both.

To overcome these limitations, we adopted a single-cell RNA sequencing (RNA-seq) approach to investigate the dynamics and heterogeneity of the stress response of stem and progenitor cells to IFNα treatment. Bone marrow cells were collected from IFNα-treated mice 3, 24, and 72 h after treatment. Cells from PBS-treated mice were included as a control. LK (Lin⁻ and cKit⁺) cells were sorted to capture a wide spectrum of the HSPC transcriptional landscape. Because HSCs are much less frequent than the other populations in the LK gate, we enriched the sorted LK samples with LK CD150⁺ CD48⁻ CD34⁻ cells at a fixed ratio to the number of LK cells to guarantee sufficient numbers of stem cells for analysis (Fig S1A). Inter-animal heterogeneity of the inflammatory stress response was addressed by performing cell hashing, for which cells from each biological replicate and time point were labeled with a unique hashtag antibody (Fig 1H and the Materials and Methods section). The cells from all four experimental time points (control, 3, 24, and 72 h) were sequenced simultaneously. The resulting dataset was processed using a computational pipeline that we designed specifically for poststimulation single-cell time series (Fig 1I). First, the cells that did not meet the quality control standards (e.g., doublets, dying cells) were removed, resulting in a total count of 1,600–3,500 cells per time point (see the Materials and Methods section). Next, clustering of the cells identified 14 different clusters in the control subset (Fig S1B). Using known marker genes for HSPC types, a scoring for stemness and a comparison with a previously published dataset (Nestorowa et al, 2016), eight different cell types could be distinguished in the control subset (Figs 1J and K and S1C–E). To compare the identified cell with to their equivalents in the treatment subsets, the cell type labels were transferred from the control subset to the other three treatment time points (Fig 1I and the Materials and Methods section). To do this without the treatment affecting cell type labels, data from all four time points were “batch-corrected.” In the batch-corrected datasets, cells in posttreatment timepoints were labeled based on their nearest neighbors in the control dataset. In addition, the batch-corrected data were used to produce a UMAP where cells are located based on their cell type identity, rather than their experimental timepoint (see difference between the two bottom UMAPs in Fig 1I). It is important to point out that the batch-corrected data were only used for label transfer and visualization purposes, but for all downstream analyses, the original data were used. Marker gene expression confirmed the cell type labels in the response time points (Fig S1F). Analysis of the hashtags showed that the biological replicates in each time point had comparable abundances of each cell type label (Fig S1G). In addition, expression analysis of the IFN α/β receptor (Ifnar1 and Ifnar2) confirmed that all clusters expressed the IFN α/β receptor and thus were able to directly respond to IFNα (Fig S1H–K). As a first measure of the inflammatory response, the expression of ISGs was scored (Fig 1L and M). At 3 h, all clusters underwent a change in their ISG expression, indicating that the whole HSPC compartment sensed the IFNα treatment (Fig 1L). However, the data also indicate great heterogeneity in the ISG response between and within clusters with HSC clusters showing the biggest change (Fig 1M).

Inflammation response is defined by global and cluster-specific changes in gene expression

To gain a comprehensive understanding of all genes that characterize the IFNα response, we next performed differential gene expression analysis. Differentially expressed genes (DEGs) were selected between the control subset and every treatment time point individually to get a set of response genes from any stage of the response (see the Materials and Methods section). The analysis identified a total of 2,501 significant response genes. Expression profiles of the response genes showed that the expression of some genes changed globally (e.g., Cox7c), whereas the expression of others was more specifically changed in few or one cluster (e.g., Sec61g and Mnda) (Fig 2A). To investigate in which cluster(s) genes were changing the most, the top 500 most significant response genes were scored for the total expression change in each cluster (change score; see the Materials and Methods section). After calculating the total change for each cluster, the response genes were categorized into 14 different groups using hierarchical clustering (Fig 2B), confirming a wide variety of globally responding genes (groups 1–5), and cluster-specific responding genes (groups 6–14). To exclude that the differences between clusters were a result of completely different expression profiles, the similarity between a cluster’s expression profile and the expression profile in the whole dataset was calculated (Fig S2A). The results indicate that most expression profiles follow similar patterns in all clusters.

Gene ontology (GO) enrichment analysis of the global response gene signatures in groups 1–4 revealed an overrepresentation of terms associated with translation and metabolism (Figs 2C–E and S2B–D). This is in line with reports of HSPCs undergoing massive changes in the metabolism under inflammatory stress (Karigane & Takubo, 2017). In addition, global response genes from group 5 (Fig S2E) showed enrichment for categories associated with immune response and response to type-I IFN, further supporting the ISG expression data (Fig 1K), indicating that all cells sense the changes in IFNα levels. The GO enrichment terms for other gene groups are shown in Fig S2F–L. Interestingly, expression changes in the HSC-enriched groups 12 and 14 were also associated with immune response and response to type-I IFN (Figs 2E and S2K). However, these changes were different from the changes in group 5, suggesting an HSC-specific immune response, which is different from the immune response in progenitors. HSC-enriched groups 12 and 14 also included GO terms such as regulation of T cell activation and antigen processing and presentation (Figs 2E and S2K), which correspond to the newly identified role of HSCs as immunomodulators (Hernández-Malmierca et al, 2022), which would be strengthened under inflammation. Besides global and HSC-specific response genes related to immune response, change score analysis also identified groups of response genes enriched in committed progenitors related to progenitor-specific processes. For example, erythroid and basophil/mast cell/eosinophil progenitor-enriched groups 9 and 11 showed an overrepresentation of processes related to erythrocyte differentiation terms and myeloid development (Fig S2H and J). Change scores in groups 8 and 10 were largest for myeloid progenitors and connected with biological processes such as phagocytosis, myeloid leukocyte-mediated immunity, and stem cell differentiation, which are characteristic functions of this cell type (Figs 2D and S2I). Thus, with the change score analysis both global and cluster-specific signatures were identified. The analysis highlights that HSCs are the major responders to inflammation in the HSPC compartment and both global and HSC-specific inflammation signatures are present, indicating heterogeneity in the inflammatory response between the clusters.

The pseudotemporal ordering of cells enhances the resolution of gene dynamics

The change score analysis gave an overview of all response genes without taking into account the detailed dynamics of the response. Therefore, in the next step, the expression dynamics of the response genes were explored. When zooming into the expression of individual genes in time, different expression patterns were observed for genes within the same group. For example, the responses of Aldh1b1 and F13a1 were both assigned to group 8. However, the temporal expression dynamics differed between those genes; whereas Aldh1b1 was up-regulated with a peak at 3 h and a full recovery to original (control) expression levels at 24 h (Fig 2F), F13a1 expression steadily went down (Fig 2G). To improve the characterization of the expression patterns, we wanted to leverage the single-cell resolution of our dataset. Therefore, we aimed to construct a pseudotime axis in the gene expression space to describe the inflammatory response. In datasets covering a developmental process or disease progression, the asynchrony of cells can be leveraged to infer a pseudotemporal ordering of cells (using methods such as diffusion pseudotime [Haghverdi et al, 2016], Monocle [Qiu et al, 2017], etc.). These methods are generally based on cell neighborhood relations, with the assumption that the further apart (in Euclidean, diffusion space, etc.) two cell states are, the longer the typical transition time between them, hence longer pseudotime (Fig 3A). However, this assumption is violated for the type of post-drug treatment time series data we have here, where the largest transcriptional change is observed shortly after stimulation but (presumably) diminishes as cells relax to a more control-like state over a longer time (Fig 3A). In our datasets, cells from different experimental time points generally appear either completely intermingled or completely disconnected from one another when viewed over the first principal components (Fig S3A). This renders the unsupervised, distance-based pseudotime methods unsuitable for detecting the temporal order of response dynamics. Thus, to get a higher temporal resolution of the response progression from the four time points, we used a semi-supervised (using the experimental time labels) approach that finds a pseudotime axis that best correlates with the cells experimental time labels. We find a transformation of the cells’ expression matrix that reconstructs the experimental time labels plus an error term that we try to minimize (Fig S3B and the Materials and Methods section). In essence, this approach assigns a positive/negative weight to each gene (depending on its correlation or anticorrelation_across all cell types_ with the experimental time labels), in such a way that the resulting weighted sum of the gene expression profile of a cell approximates its experimental time labels up to an error term, which we interpret as the asynchrony of treatment response among different cells captured at the same experimental time point. Using this approach implies that there are at least a few genes in the data, a linear combination of which can be used to differentiate the different time points. Interestingly, the genes given highest positive weights by the model turn out to be often among the global changing ones, and associated with the immune response, whereas the genes assigned negative weights display enrichment in processes such as translation (Fig S3C–G). In Fig S3H, we see the non-smoothed counterpart, thus sparser expression matrix of the same genes as that of Fig S3C, which confirms the robustness of the inferred patterns. We further confirm the agreement of actual time series dynamics with that of the reconstructed pseudotime pattern for a few show-case genes (Fig S4A–G). In addition, we introduce a simulation dataset on which we demonstrate the properties of the inferred pseudotime order, and use it for evaluation of the method’s performance (the Materials and Methods section and Fig S5A–M).

Using the response pseudotemporal order of the cells, we moved from the four time points’ discreet view of the data to a continuous axis (Fig 3B) over which the different expression dynamic patterns that follow IFN⍺ treatment were explored. Response genes were categorized into 16 patterns based on their pseudotemporal dynamics using hierarchical clustering (Figs 3C and S3H). Each pattern represents a cluster of genes with similar expression dynamics after IFN⍺ injection, which can be roughly subdivided in up-regulation (patterns 1–9 and 16) and down-regulation (pattern 10–15) patterns. The GO enrichments for the gene patterns are shown in Fig 3D–G. The patterns showed a broad diversity in the speed of sensing, response, and recovery. Whereas most patterns reach a steady-state plateau towards the end of the pseudotime axis, a few (patterns 5, 8, and 12) do not, and imply an ongoing trend of changes beyond the covered 72 h posttreatment genes in patterns 6, 7, and 9 show quick sensing, response, and recovery associated with immune, IFN, and viral processes (Fig 3F). Other patterns resembled a similar fast increase but with slower recovery, as seen for pattern 3, which is enriched for metabolic processes (Fig 3E). In addition, the heatmap in Fig 3C showcases a variety of gene dynamics that were different from the rapid response and recovery IFN-response (patterns 6, 7, and 9), such as a sustained up-regulation in pattern 1, which was associated with translation and other biosynthetic processes (Fig 3D). In contrast, several gene patterns encompassed genes that were down-regulated (gene patterns 10–14). After the initial decrease in expression, many of these genes failed to recover to initial expression levels. Most of these genes were linked to myeloid development and differentiation (Fig 3G). This would suggest alterations in myeloid differentiation upon the IFN⍺ treatment, an observation described for many other pro-inflammatory cytokines (Matatall et al, 2014; Pietras et al, 2016; Yamashita & Passegué, 2019) but IFN⍺.

Response pseudotime reveals a landscape of gene dynamics in HSPCs after IFN⍺ treatment

To decipher the dynamic changes in the inflammatory response in the different clusters, we combined the information on how response genes changed their expression (Fig 3C) with whether these changes were global or cluster-specific (Fig 2B). The result condenses the plenitude of information in the complete single-cell time series into a single visualization, which considerably eases the search for (groups of) biologically relevant genes (Fig 4A).

The most commonly found patterns among the groups are patterns 1, 2, 7, and 9, showing a fast up-regulation combined with a partial (patterns 1 and 2) or full (patterns 7 and 9) recovery. Genes from these patterns are mainly related to immune responses, highlighting that the sensing and response to IFN⍺ is fast and present in all groups and clusters. The full recovery in the most common patterns 7 and 9 revealed by our single-cell experiment indicates that the general immune response does fully recover within 72 h, a fact that bulk experiments had not captured before.

Interestingly, other patterns of fast sensing and response followed by sustained up-regulation (patterns 1 and 3) were enriched in the groups with a global signature (groups 1, 2, 3, and 4). Several of these genes were ribosomal (e.g., Rpl35 and Rps27; Fig 4B), suggesting that biosynthetic activity increased in most HSPCs early in the treatment and remained active even in the recovering phase, when proliferation levels are back to homeostasis again (72 h; Fig 1B). In addition, several genes after the sustained up-regulation pattern were metabolic. Oxidative phosphorylation (OXPHOS) genes and mitochondrial enzymes (e.g., Atp5e and Cox7c; Fig 4C) showed these patterns of prolonged up-regulation. On the other hand, glycolytic genes Hk2 and Pgk1 showed a quick response and recovery (Fig 4D). Thus, in contrast to previous reports suggesting a binary (on/off) switch between glycolysis and OXPHOS (Suda et al, 2011), our data suggest that an initial up-regulation of glycolytic and a sustained up-regulation of OXPHOS genes go hand in hand in inflammation responding HSPCs.

In contrast to the heterogeneity in dynamics observed globally, HSC-enriched groups (groups 5, 12, and 14) are mainly enriched with gene patterns that increase very early after treatment and quickly return to homeostatic levels (patterns 6, 7, and 9). Examples of such genes are Irf7 and Irf9 in group 5 (Fig 4E) or Ifit2 and Sca-1 in groups 12 and 14 (Fig 4F and G). Thus, most of the HSC-enriched groups follow rapid sensing, responding, and recovery dynamics, with most of the gene changes preceding the peak in proliferation response in these cells. In addition, most of the HSC-enriched dynamic changes are within gene groups linked to IFN and immune response, again highlighting the specific, fast HSC-specific immune response.

In contrast to HSC-specific groups, committed progenitor-specific groups (8, 9, 10, and 11) were strongly enriched in genes that exhibited persistent down-regulation (patterns 10, 11, 12, and 14). In the myeloid progenitor-specific groups 8 and 10, many of these genes were associated with myeloid cell differentiation and functional programs, e.g., Csf1r; Irf8 (Fig 4H and I), suggesting reduced myeloid differentiation in the myeloid progenitor clusters.

In summary, response pseudotime has shed light on the heterogeneity in gene dynamics in the HSPC compartment during the induction of inflammation. Whereas global groups encompass diversity in gene patterns, cluster-enriched gene groups show far less variation and more specificity, with HSCs being the fast responders and recovers, whereas committed myeloid progenitors showing sustained down-regulation of genes.

Single-cell abundance analysis shows myeloid depletion and HSC enrichment after IFNα treatment

To investigate whether reduced transcriptional programs for myeloid differentiation and function upon IFN⍺ treatment also impacted the size of the progenitor compartment, we performed a differential abundance analysis on the level of clusters and neighborhoods of cells. We applied the Milo algorithm (Dann et al, 2022) which models cellular states as overlapping neighborhoods on a KNN graph rather than relying on clustering cells into discrete groups (see the Materials and Methods section). At a false discovery rate of 10%, we could observe multiple neighborhoods that were differentially abundant (Fig 5A). Neighborhoods received a cell-type label based on the most predominant cluster in the neighborhood. Even though most progenitor-enriched clusters showed a reduction at 3 h, most returned to normal by 24 or 72 h, except for the most differentiated myeloid progenitors (Myel. prog. #2 and Myel. prog. #3), which were significantly reduced, even at 72 h posttreatment (Figs 5B and S6). The abundance of HSCs increased in HSCs #2 posttreatment, most significantly at 3 h (Figs 5B and S6). This can be explained by our strategy of FACS enrichment of rare cell types (the Materials and Methods section) and the fact that the total number of cells captured at the 3-h time point was smaller (the Materials and Methods section) compared with the other time points, giving the enriched HSC population a higher abundance in the 3-h sample. Our abundance analysis in the HSPC compartment also showed that acute IFN⍺ treatment resulted in a sustained significant reduction in the most committed myeloid progenitors over the whole time course of the response. This is in contrast to the current notion in the field claiming that the decreased frequency of LS⁻K (comprising myeloid, erythroid, and megakaryocytic progenitors) and concurrent increase in LSKs (comprising HSCs and LMPPs) upon IFNα stimulation (Fig 1B–E) is mainly the result of contaminating myeloid progenitors that have reacquired Sca-1 expression (and would fall into the LSK gate) (Pietras et al, 2014; Kanayama et al, 2020). However, when analyzing Sca-1 gene expression in our dataset, the absolute gene expression of Sca-1 in the myeloid progenitors never exceeds the gene expression levels of the HSCs, even though there is a relative change in each of the clusters (Fig 5C and D). Thus, this unbiased (i.e., transfer of cell type labels from the control, based on the expression of several genes) investigation of the different clusters and their abundances identifies a true change in myeloid population size and not a shift in populations because of marker change in response to inflammation.

Myeloid depletion coincides with changes in transcriptional programs

The reduction in the abundance of the myeloid progenitors could be caused by an increased egress of these cells from the bone marrow, a loss of these cells because of cell death, or a reduction in differentiation towards myeloid progenitors. Myeloid progenitors were not observed in the blood at any time point after IFNα treatment (Fig S7A). However, the number of myeloid progenitors in the spleen decreased at 24 h (Fig S7A), in line with the reduction observed in the bone marrow (Fig 5A and B). This suggests that the reduced abundance of the myeloid progenitors in the bone marrow is not because of increased egress into the blood or spleen. To investigate whether reduced levels of myeloid progenitors in the bone marrow were the result of increased cell death, gene patterns of pro-survival genes (Bcl2, Birc2, and Birc5) were analyzed and found to be decreased after IFNα treatment (Fig 6A). Yet, in Bax^−/−Bak^−/− double-knockout mice, in which cells are unable to undergo apoptosis because of the loss of the pro-apoptotic proteins Bax and Bak, a similar reduction of myeloid progenitors was observed in the bone marrow as in wild-type mice (Fig 6B), indicating that apoptosis was not the reason for the reduction in myeloid progenitors. This does not exclude the involvement of other forms of cell death, like necroptosis and pyroptosis. Expression of necroptosis-related genes (Fig 6C) was not significantly altered in any of the cell types; however, expression of some pyroptosis-related genes (Fig 6D) such as Casp1 and Casp4 (Figs 6E and S7B–D), was increased in the myeloid progenitor subpopulations. Together, these data suggest an early increase in pyroptosis rate in the myeloid progenitors compared with the other cell types in response to the treatment, possibly resulting in a reduced abundance of these cells.

To investigate whether insufficient production of new cells might also play a role in the sustained reduction in myeloid progenitors, the expression of genes involved in myeloid lineage priming was analyzed. By calculating a gene score for known myeloid transcription factors that have been reported to impact both stem and progenitor cells, a global down-regulation in the myeloid program (including priming in HSCs and the more differentiated LMPPs) was present in the early stages of the IFNα response (Fig 6F). In addition, a reduction in cell cycle and purine nucleotide synthesis genes was observed in myeloid progenitors, suggesting that both myeloid differentiation in HSCs and LMPPs and cell production in myeloid progenitors, could be affected (Fig S7E–H).

Along with the reduction in myeloid differentiation related genes, neutrophil, and monocyte transcriptional signatures were also down-regulated in myeloid progenitors over the pseudotime axis (neutrophil: myel. prog. #1; monocyte: myel. prog. #3) (Fig 6G and H), suggesting continuously reduced differentiation of myeloid progenitors to mature myeloid cells. This was confirmed by a gradual decrease in the number of neutrophils in the blood over time (Fig 6I). The number of monocytes in the blood also decreased in response to treatment, but already recovered at 72 h (Fig 6J). Similar trends were observed in bone marrow and spleen; however, no significant changes were detected (Fig S7I–L). Taken together, these in vivo cell analyses and gene expression data indicate an increase in pyroptosis signature combined with reduced myeloid differentiation both at the stem cell level, and in committed progenitors, accompanied with a reduction in the cell production machinery in progenitors, resulting in reduced myeloid progenitors in the bone marrow, and lower levels of mature myeloid cells in the blood (Fig 6K).

Limitations of the study

By creating this resource, we have established a foundation for exploring the functional properties of the molecular signatures of HSPCs under inflammatory stress. Whereas correlations inferred from expression data offer hypotheses for regulatory mechanisms, experimental testing is crucial for confirming these hypotheses. However, inflammation-induced marker changes on HSPCs are still a limitation to perform cell cluster-specific functional follow-up studies. Recently developed CITEseq analysis would allow us to identify better markers to distinguish the different stem cell and progenitor populations under inflammation to allow better isolation of these populations for future functional analysis. Moreover, additional layers of regulation at the chromatin and protein level impact and control cell behavior, which may manifest as unappreciated heterogeneity and dynamic properties. Integrating our time series data with methylation, chromatin accessibility or proteome analysis will achieve further comprehensive modeling of the regulatory and signaling networks coordinating the stress response.

We acknowledge the possibility of changes in cell abundance as a result of our rare HSC population enrichment. However, we emphasize that the significant decrease in myeloid progenitors persisted consistently across all time points, indicating a true biological change rather than an artifact (Figs 5A and B and S6). Whereas we recognize the potential influence of our enrichment for HSCs on cell abundance of other populations, the abundance reduction was specific to the myeloid population rather than a uniform reduction across all populations (as would be expected in compensation for HSC numbers). We further validated the reduction of myeloid progenitors in both blood and spleen (Fig S7); this further confirms our HSPC single-cell abundance results.

Our combined abundance and response pseudotime analysis suggest a link between myeloid progenitor reduction in the bone marrow and the HSCs reduced myeloid differentiation bias and a potential enhancement of cell death mechanisms in the myeloid progenitors, relative to the other HSPC populations.

Validation of these findings at the transcriptional level is however not straightforward because of the multitude of posttranscriptional regulatory mechanisms involved until a specific biological function (e.g., enhanced cell death) is manifested in the system. We also sought to check if the post-inflammation decline in differentiation bias towards myeloid progenitors could be confirmed by a comparison of RNA-velocity of the cells at the different time points. We estimated RNA-velocities of the cells in the control data set (Marot-Lassauzaie et al, 2022). However, meaningful comparison with the cell velocities of later time point datasets was hindered by the yet insufficient robustness of state-of-the-art velocity parameter inference and the poor quantification of unspliced versus spliced mRNA, especially when using the 3′-end gene expression short reads captured by the 10X sequencing platform (Marot-Lassauzaie et al, 2022). Emergence of more reliable approaches for estimating cell velocities could be a valuable asset as an orthogonal dynamical analysis approach to time series and pseudotime analysis, before we (or other researchers) turn to perform further costly and challenging functional and experimental validations in the future.

Mouse models

All animal experiments were approved by the local Animal Care and Use Committees of the German Regierungspräsidium Karlsruhe für Tierschutz und Arzneimittelüberwachung. Mice were kept under specific pathogen-free conditions in ventilated cages (ICV) in the animal facility of the German Cancer Research Center (DKFZ). Mice used for experiments were between 10–20 wk old at the beginning of the respective experiments. SclCreERT bax^−/−bak^−/− mice were on a C57BI/6 background (Takeuchi et al, 2005), and treated for 5 d with 2 mg/d tamoxifen. IFNα treatment was started 4 wk post tamoxifen treatment. C57Bl/6 (WT) mice were bred at the DKFZ animal facility or bought from JANIVER lab. Mice were euthanized by cervical dislocation according to German guidelines.

IFNα treatment of mice

Mice were injected subcutaneously with 50,000 international units (IU) of recombinant mouse IFNα per 20g mouse (Milteny Biotech). Recombinant mouse IFNα was diluted in PBS and control mice were injected with 100 μl PBS.

Isolation of BM, spleen, and blood for flow cytometry analysis

Blood was collected from the vena facialis by sub-mandibular bleeding into EDTA-coated collection tubes. Blood was either analyzed automatically with a Hemavet cell counter (Drew Scientific) or stained for flow cytometry after initial RBC lysis by incubation with ACK lysis buffer for 20 min. Cells were stained for Ter119, CD4, CD8, CD11b, CD11c, Ly6G, B220, F4/80. BM cells were isolated from the femur, tibia, hip bone, and spine by bone-crushing. Splenocytes single-cell suspension was obtained by mashing the spleen through a 40 μm EASYstrainer (greiner bio-one). After ACK, lysis BM cells and splenocytes were stained using antibodies for CD117 (cKit), Sca-1, CD150, CD48, CD34, CD16/32, and lineage antibodies (CD4, CD8, CD11b, Gr-1, B220, and Ter119). For the BrdU incorporation assay, BrdU (18 mg/kg; Sigma-Aldrich) was administered i.p. for 14 h before harvesting the BM. The BD Pharmingen BrdU Flow Kit protocol was used to stain for BrdU. For flow cytometry analysis the LSR Fortessa or LSRII were used (BD Biosciences). Flow data were analyzed using BD FACS DIVA v8.0.1 and Flowjo (v10).

FACS sorting

For FACS sorting of single cells, BM cells were isolated and RBC lysed as described above. This was followed by lineage depletion using a lineage antibody cocktail against CD4, CD8, CD11b, B220, Gr-1, and Ter119 and incubation with Dynabeads Magnetic Beads (Invitrogen). Lineage-depleted BM cells were stained with Zombie Yellow viability dye (BioLegend) followed by incubation with the following antibodies: CD117, Sca1, CD150, CD48, CD34, and lineage antibodies (B220, CD4, CD8, Ter119, Ly6G, CD11b, CD11c, and F4/80) together with one of the hash antibodies (TotalSeq-A0301 anti-mouse Hashtag 1 Antibody, TotalSeq-A0302 anti-mouse Hashtag 2 Antibody, TotalSeq-A0303 anti-mouse Hashtag 3 Antibody, TotalSeq-A0304 anti-mouse Hashtag 4 Antibody) (TotalseqA antibodies; BioLegend). The four biological replicates of each time point were stained with one of the four unique hash antibodies. Cells were sorted using a FACSAria Fusion or FACSAria II equipped with a 100 μm nozzle (BD Biosciences).

Single-cell RNA library preparation and sequencing

HSPC single-cell RNA-seq was performed using the 10X Genomics platform. The Chromium Next GEM single cell 3′ reagent kits v3.1 were implemented to prepare the libraries, following the official instruction manual (https://www.10xgenomics.com/support/single-cell-gene-expression/documentation/steps/library-prep/chromium-single-cell-3-reagent-kits-user-guide-v-3-1-chemistry). Briefly, 10,000 Lin⁻ cKit⁺ cells were sorted and enriched for HSCs by sorting additional 3,000–4,000 Lin⁻ cKit⁺ CD150⁺ CD48⁻ CD34⁻ cells. Cells were super-loaded according to the manufacturer’s instructions up until the cDNA amplification step. 1 μl/sample of HTO primers was spiked into the cDNA amplification PCR, and cDNA was amplified according to the 10x Single Cell 3′ v3.1 protocol aiming for a targeted cell recovery of 500–6,000 cells. After PCR, cDNA cleanup was performed by using SPRI to separate the HTO-derived cDNAs (in the supernatant) from the mRNA-derived cDNAs (retained on beads). The cDNA fraction was processed according to the manufacturer’s protocol to generate the transcriptome library. The quality of the obtained cDNA library upon adapter ligation and sample index PCR was assessed on an Agilent Bioanalyzer High sensitivity chip. Library sequencing was performed on the Novaseq 6000 Illumina sequencing platform.

Filtering longitudinal single-cell RNAseq dataset

The cellranger pipeline (version 3.1.0) was used to align all reads to the mm10 genome and count the coverage of each gene in each cell. Based on the hashtag barcodes, cells were assigned to their corresponding time point (control, 3, 24, or 72 h) and batch (four batches per time point). Cells with multiple barcodes (multiplets) or missing barcodes (negatives) were removed from the dataset. In the resulting count matrix (cells × genes), cells with a high amount of mitochondrial genes (>5%) or a low amount of unique genes (<700) were filtered out. After the filtering steps, the following number of cells was present in each of the respective time points: control—2,474, 3 h—1,661, 24 h—3,462, 72 h—2,449 (10,046 cells in total).

Clustering and cell type annotation

The 500 most highly variable genes (HVGs) were identified in the control subset using analytic Pearson residuals (Lause et al, 2021). The control subset was subsetted for the 500 HVGs, and the counts were L2 normalized. Next, a neighborhood graph was computed using 10 out of 50 principal components and the 15 nearest neighbors. The Leiden algorithm (resolution = 0.8) identified 14 distinct clusters in the control subset (Traag et al, 2019). Each cluster was appointed to a cell type based on (1) DEGs between the cluster of interest and all other clusters, (2) the expression profiles of the HVGs, (3) known marker genes, and (4) correlation with cell types in a previously published dataset of the HSPCs (Nestorowa et al, 2016).

Label transfer and UMAP representation

We identified the top 2,000 HVGs in each subset and subsetted the complete dataset with the combined list of HVGs. Afterward, the dataset was L2 normalized, and the different subsets were integrated using Scanorama (Hie et al, 2019). All 100 Scanorama-reduced dimensions were used to calculate a neighborhood graph (nearest neighbors = 15). A two-dimensional UMAP representation was computed using the neighborhood graph. To transfer the cell-type labels from the control subset to the response subsets (3, 24, and 72 h), cells in the response subsets would adopt the cell type label that was most common among their 15 nearest neighbors (Euclidean distance) in the control subset. The integrated data were only used for label transfer and visualization purposes. For other downstream analyses, we use the filtered-only dataset. In this dataset, we removed the Ba/MC/eos. cells and monocytes because of the small number of cells assigned to those cell types (10 and 52, respectively).

Calculating gene set scores

The filtered dataset was L2 normalized and scaled to unit variance and zero mean. The ISG score was calculated by subtracting the average expression of a random set of reference genes from the average expression of about 400 known ISGs (Scanpy function score_genes). Similarly, the stemness (Giladi et al, 2018), necroptosis (GO:0070266), pyroptosis (GO:0070269), myeloid TF (Kwok et al, 2020), monocyte and neutrophil differentiation, cell cycle (Giladi et al, 2018), and purine nucleotide synthesis (Vogel et al, 2019) score were calculated. Genes for each signature are available as a .csv file. Necroptosis and pyroptosis gene sets were retrieved from the Mouse Genome Database, Mouse Genome Informatics, The Jackson Laboratory. World Wide Web (URL: http://www.informatics.jax.org). (The data were retrieved in the year 2022).

Differential abundance analysis

We used the R package Milo to perform an abundance analysis on the L2 normalized, filtered dataset (Dann et al, 2022). A neighborhood graph was built using 30 out of 100 of the Scanorama-reduced dimensions (see Label transfer and UMAP representation) and 30 nearest neighbors. Afterward, we followed the steps described in the accompanying tutorial (Milo example on mouse gastrulation dataset) for each response subset (3, 24, and 72). In each analysis, the control subset served as the reference, to which the response subset would be compared.

Identifying response genes

We used the edgeR-LRT method in the Libra R package to find the DEGs between the control and any of the response subsets, in each cluster. We considered only DEGs with an adjusted P-value higher than 0.05 and a log-fold change higher than 1 in at least one cluster. For the downstream analyses of the response genes, we consider only the 500 DEGs with the highest significance. In case a DEG is found in more than one cluster and/or time point, we take only the lowest P-value into consideration.

Change score

We L2 normalized the filtered dataset per cell. For each response gene (i) we took the mean expression (μ) in each cluster (j) per time point (t). The expression change was calculated as the absolutes sum of the derivative of the mean expression across all time points (here m = 3 for control-3, 3–24, 24–72 h). Thus the change score (c_i,j) per cluster for each response gene:ci,j=∑t=1m|μi,j,t+1−μi,j,t |(1)

The result is a matrix with change scores per cluster for each of the response genes. We applied hierarchical clustering and grouped the response genes into 14 groups by setting a threshold at the cophenetic distance of 3 (Scipy function cluster.hierarchy.linkage and cluster.hierarchy.fcluster).

Similarity score

We define the similarity score between the cluster-specific expression profiles of gene i and the expression profile of the same gene in the complete dataset in two steps. First, both the expression of the cluster and the complete dataset were scaled between 0 and 1 by min–max normalization (with n = 4 indicating the number of time points in the time series).x_i,j,t=μi,j,t − mint μi,j,tmaxt μi,j,t − mint μi,j,t(2)

Note that max_t μ_i,j,t and min_t μ_i,j,t indicate the max/min of the mean expression for gene i cluster j among all time points t. Second, we calculate a dissimilarity score between the cluster and the whole data set for each response gene by subtracting the cluster-specific expression changes from the expression changes in the complete dataset. Third, the magnitude of these differences were summed and normalized by the number of time points. Finally, we turn the dissimilarity into a similarity score (s_i,j) by subtracting it from 1 such that 1 presents complete similarity to the average behavior of all cells in the dataset.si,j=1 − ∑t=1n|x_i,j=all,t − x_i,j,t|n(3)

The result is a matrix with similarity scores per cluster for each of the response genes.

Pseudotemporal ordering of cells during response

We opt to find a pseudotime axis that correlates with the actual arrow of time. Thus, we look for a transformation (W of size [G,1]) of the expression data from all time points that reconstructs the experimental time point of each cell with minimal error (ε):X∗W=T+ε(4)

Here, X (of size [N,G]) is the filtered count matrix after per-cell L2 normalization, subsetting for all response genes and per-gene standardization (i.e., setting the mean to zero and variance to 1). The per cell L2 normalization is part of our single-cell data preprocessing routine (also see Haghverdi et al [2018]) on our preference of L2 normalization rather than size factor L1 normalization). The per-gene standardization is for correct calculation of the X^T X matrix (also necessary for other methods such as principal component analysis), such that the different features become comparable; for example, genes with overall high expression do not artificially seem highly correlated to other genes. Thus, by data standardization, we ensure mean- and variance-independent calculation of the weights in the W matrix. T is a vector (of size [N, 1]) with an (experimental) time assignment for each cell, created by taking the experimental time points (control, 3, 24, and 72 h) and converting those to 0, 1, 2, and 3 values respectively (as an alternative to such equal importance separation of the time points, one could consider using the actual time values on a log-scale). By this choice, we imply equal interest in separating each time point from the subsequent time point, whereas e.g., using the actual time values (0, 3, 24, and 72) would mean a larger emphasis on finding separation between the last two time points with the largest distance (i.e., 72−24 = 48).

In this work we have used the regression model (without an intercept) on standardized expression values and time scale choice as described above. Nevertheless, one could also consider other equally well justifiable variations, for example, a regression model including an intercept term. Note that the notion of pseudotime, which is generally based on a transformation of cell states’ profiles (e.g., gene expression), provides only a flexible estimation of the dynamical progression order of the cells and is not an absolute or unique quantity (e.g., when using slightly different transformations) as discussed in supplementary note 6 of Haghverdi et al (2016), and also demonstrated on simulation data (Fig S5J–M). The least squares analytical solution for W (which minimizes ε^T ∗ ε) is given by:W=(XTX)−1∗XT∗T(5)

We used the expression matrix X with the size of 10,046 cells and 2,501 genes and the cells’ corresponding time labels to solve the above linear regression problem. We note that, to avoid over-parametrization and to ensure the identifiability of the solution, the number of cells has to be larger than the number of genes. After W has been retrieved, a pseudotime coordinate can be calculated for each cell by the following:PT=X ∗ W(6)Where PT (of size [N, 1]) is the vector with the pseudotime coordinate for each cell. For further downstream analyses, the cells were ordered based on their pseudotemporal assignment. Note, the genes which the model forms a linear combining of, do not necessarily need to be linear or monotonic themselves. In fact, on simulation data, we show the model is able to reconstruct the pseudotime, (up to an acceptable correlation with the true time orders), even in a case where all genes are non-monotonic (complex) in the data (Fig S5I).

To cluster the response genes based on their pseudotemporal expression pattern, we smoothed the expression of these genes by taking the mean expression per 100 cells over the pseudotime axis. The expression values were scaled between 0 and 1 by min–max normalization. The 500 response genes were then clustered into 16 pseudotemporal expression patterns, using hierarchical clustering (threshold at cophenetic distance of 5.2).

Simulation of time series

For the discrete time simulation, we sampled 1,000 cells from four different discrete timepoints: 0, 5, 10, and 20 (250 cells per timepoint). To model the asynchronous progression of time, we added a Gaussian noise (SD 1.5) to the discrete sampling time of each cell, giving us each cell’s hidden pseudotime t. For the continuous time simulation, we uniformly sampled cells pseudotime t between 0 and 20. The expression dynamics of 20 different genes were then simulated using either a single sigmoid function: S(t|α,β)=11+e((−t−α)/β) with randomly sampled offset α and width β (monotonic genes) or a sigmoid function deactivated at a switching time point t_switch : S(t|α,β) - S(t - t_switch|α,β) (complex genes). The genes are then min–max normalized and Gaussian noise SD of 0.05 is added to each cell’s expression to simulate biological variance (see also Fig S5).

Gene expression over pseudotime

Expression profiles of individual genes in response pseudotime (such as Fig 4B) were derived using a combination of bin smoothing and bootstrapping. To find the expression profile in the complete dataset, bin smoothing with a 600-cell window size was performed on a sample of 50% of the cells in the dataset. This was repeated 20 times to find the mean expression, which defines the expression profile. The 95% confidence intervals were calculated by multiplying the standard error with 1.96 and subtracting or adding to the mean. For the cluster-specific expression profiles of individual genes, a 50-cell window size was chosen instead, because of the smaller number of cells in each cluster.

Gene score in pseudotime

The gene set score profile in response pseudotime was calculated using a combination of locally weighted least squares regression (LOESS) smoothing and bootstrapping. For each cluster, LOESS smoothing with a first-order regression model was applied to 50% of the cells. This was repeated 30 times. The score profile was derived by taking the mean and 1.96 times the standard error for the 95% confidence intervals.