Neutrophils and emergency granulopoiesis drive immune suppression and an extreme response endotype during sepsis

Sepsis arises from diverse and incompletely understood dysregulated host response processes following infection that leads to life-threatening organ dysfunction. Here we showed that neutrophils and emergency granulopoiesis drove a maladaptive response during sepsis. We generated a whole-blood single-cell multiomic atlas (272,993 cells, n = 39 individuals) of the sepsis immune response that identified populations of immunosuppressive mature and immature neutrophils. In co-culture, CD66b+ sepsis neutrophils inhibited proliferation and activation of CD4+ T cells. Single-cell multiomic mapping of circulating hematopoietic stem and progenitor cells (HSPCs) (29,366 cells, n = 27) indicated altered granulopoiesis in patients with sepsis. These features were enriched in a patient subset with poor outcome and a specific sepsis response signature that displayed higher frequencies of IL1R2+ immature neutrophils, epigenetic and transcriptomic signatures of emergency granulopoiesis in HSPCs and STAT3-mediated gene regulation across different infectious etiologies and syndromes. Our findings offer potential therapeutic targets and opportunities for stratified medicine in severe infection. Knight and colleagues report altered granulopoiesis and increased frequency of immature neutrophil subsets with immunosuppressive properties in a subset of patients with sepsis with poor outcome.

Sepsis arises from diverse and incompletely understood dysregulated host response processes following infection that leads to life-threatening organ dysfunction. Here we showed that neutrophils and emergency granulopoiesis drove a maladaptive response during sepsis. We generated a whole-blood single-cell multiomic atlas (272,993 cells, n = 39 individuals) of the sepsis immune response that identified populations of immunosuppressive mature and immature neutrophils. In co-culture, CD66b + sepsis neutrophils inhibited proliferation and activation of CD4 + T cells. Single-cell multiomic mapping of circulating hematopoietic stem and progenitor cells (HSPCs) (29,366 cells, n = 27) indicated altered granulopoiesis in patients with sepsis. These features were enriched in a patient subset with poor outcome and a specific sepsis response signature that displayed higher frequencies of IL1R2 + immature neutrophils, epigenetic and transcriptomic signatures of emergency granulopoiesis in HSPCs and STAT3-mediated gene regulation across different infectious etiologies and syndromes. Our findings offer potential therapeutic targets and opportunities for stratified medicine in severe infection.
The multiple dynamic host pathophysiological mechanisms that result in organ dysfunction following infection are incompletely understood and, while often aggregated into the clinical syndrome of sepsis, overlap with other critical illness syndromes [1][2][3] . There is an urgent need to better delineate such extreme responses to infection, given the COVID-19 pandemic and the wider, global burden of all-cause sepsis, which accounts for 11 million deaths per year and has a persistently high mortality of 20-30% 4 . Axes of immune dysregulation in sepsis have typically been examined by peripheral blood bulk transcriptomic studies, which lack the resolution to identify cell type-specific signatures [5][6][7] or single-cell interrogation of peripheral blood mononuclear cells (PBMCs) 8 , which omits neutrophils (Neu).
Data from animal models demonstrate key protective and pathogenic roles for Neu in sepsis, in some instances involving specific Article https://doi.org/10.1038/s41590-023-01490-5 Sepsis is highly heterogeneous 2 and immune suppression associated with apoptosis, epigenetic reprogramming and downregulation of activating cell surface molecules in multiple cell types 2 is a predominant feature in many patients. Sepsis subphenotypes are reported 2,[5][6][7]16 , but their relationship with mechanisms of immune dysfunction is not well defined. Sepsis response signatures (SRSs) from whole-blood transcriptomics identify dynamic response states, with assignments in the SRS1 group or a high likelihood of SRS1 group (SRSq), associated subsets 9,10 . In humans, the clinical observation of a 'left shift' in complete blood count to increased immature Neu in severe infection is well recognized 11 . Neu abundance varies between septic and non-infectious inflammation 12 , while single-cell transcriptomics indicates that they are important in the pathogenesis of severe COVID-19 (ref. 13) and acute respiratory distress syndrome (ARDS) 14 . Expansion of subsets of immunosuppressive granulocytes associate with a higher risk of nosocomial infection 15   , sterile inflammation controls after CS (n = 7) and patients with sepsis (n = 26) by joint scRNA-seq and cell-surface protein profiling; scHSPC atlas for HCs (n = 7) and patients with sepsis (n = 15) (multiome scRNA-seq and scATACseq). Validation phase: mmV in HCs (n = 11) and patients with sepsis (n = 36), including CyTOF and bulk RNA-seq of whole blood and timsTOF proteomics in plasma; WB-CID datasets (n = 1,595) with deconvolution using scWB. b, Uniform Manifold Approximation and Projection (UMAP) for scWB (HCs, n = 6; CS, n = 7; sepsis, n = 26) scRNA-seq (272,993 cells) annotated for differential cell populations. DC, dendritic cell; cDCs, classical DCs; pDCs, plasmacytoid DCs. c, UMAP for scWB scRNA-seq differential abundance in samples from patients with sepsis (n = 26) compared to HCs (n = 11), with sampled neighborhoods colored by statistical significance (spatial FDR < 0.05). Nhood, neighborhood. d,e, Beeswarm plots of differential cell abundance in scWB with cluster labels of neighborhoods depicted and compared for patients with sepsis (n = 26) versus HCs (n = 11) (d) and patients with sepsis (n = 26) versus after CS (n = 7) (e).
Article https://doi.org/10.1038/s41590-023-01490-5 with immunosuppression, differential response to steroid therapy, more severe disease and higher early mortality 5,[17][18][19] . Here we identified differences in Neu function and subset abundance during sepsis, which were a consequence of altered granulopoiesis. We demonstrated these immunosuppressive granulocytic and granulopoietic disturbances were the functional basis of the SRS1 subphenotype of sepsis.

Immature and cycling neutrophils are increased in sepsis
To generate an unbiased single-cell whole-blood (scWB) atlas of the sepsis response for all peripheral blood leukocyte populations, including Neu, we assayed freshly sampled whole blood from 26 patients with all-cause sepsis with a change in quick sequential organ failure assessment (SOFA) score of ≥2 points, indicating organ dysfunction and a physiological measure of acute illness severity score (national early warning score 2; NEWS2) of ≥7, indicating patients who required an urgent critical care response (Supplementary Table 1 and Methods); 9 sepsis convalescents (sampled 1-3 months after hospital discharge); 6 age-and sex-matched healthy controls (HCs); and 7 patients after cardiac surgery (CS) as a sterile inflammation control ( Fig. 1a and Methods).  Table 2). After clustering and annotation of major immune cell types (Extended Data Fig. 1a-c and Methods), we observed recognized hallmarks of peripheral blood sepsis immunophenotypes, including neutrophilia, lymphopenia and reduced HLA-DR expression in CD14 + CD16 − classical monocytes (cMo) 20 (Extended Data Fig. 1d).
We performed fine-resolution clustering and annotation (Extended Data Fig. 1e-i and Supplementary Table 3). RNA velocity and partition-based graph abstraction showed annotated Neu subpopulations followed the expected maturation sequence from immature MPO + Neu to PADI4 + Neu to IL1R2 + Neu to S100A8/9 + Neu to mature CD10 + CD16 + Neu (Extended Data Fig. 1j,k). We identified that degranulating CEACAM8 + Neu, S100A8/9 hi Neu, IL1R2 + Neu, PADI4 + Neu, MPO + Neu and cycling MK167 + CYP1B1 + Neu were all proportionally increased in sepsis compared to HCs, whereas all mononuclear cell subsets, except CD19 + CD38 + CD71 + plasmablasts, were reduced (Fig. 1c,d). The increase in degranulating CEACAM8 + Neu and S100A8/9 hi Neu and the reduction in all mononuclear cell subsets, was also seen in CS compared to HCs (Extended Data Fig. 2a,b), suggesting they represented nonspecific features of inflammation. By contrast, higher abundance of the immature IL1R2 + Neu, PADI4 + Neu and MPO + Neu subsets and cycling MK167 + CYP1B1 + Neu were specific to sepsis (Fig. 1c-e and Extended Data Fig. 2c,d). These findings demonstrated differential abundance of specific immature and cycling Neu subsets in sepsis.

Neutrophils in sepsis and recovery are immunosuppressive
To functionally test the immunosuppressive properties of sepsis Neu, we co-cultured bulk CD66b + Neu (isolated from fresh whole blood of patients with sepsis by immunomagnetic selection) with allogeneic CD4 + T cells isolated from healthy donor leukocyte cones at a 4:1 Neu:T cell ratio in medium supplemented with interleukin (IL)-2 and CD3/CD28 Dynabeads. After 72-96 h of culture, the fraction of proliferating CD4 + T cells (calculated using non-bead-stimulated T cells cultured without Neu as the baseline proliferative fraction of cells in each sample) and the percentage of either PD-1 + or CD69 + CD4 + T cells (calculated relative to the CD3/28 bead-stimulated, no T cells control culture, to account for donor variation) was lower when co-cultured with CD66b + Neu from patients with sepsis compared to HCs (Fig. 2a), indicating sepsis CD66b + Neu inhibited CD4 + T cell proliferation and activation. There was no difference in the percentage of live, non-apoptotic CD4 + T cells on far-red DNA staining after 72-96 h of co-culture with either HCs or sepsis CD66b + Neu (Extended Data Fig. 2e). For samples on which both single-cell sequencing and functional co-culture assays were performed, we correlated MPO + Neu, PADI4 + Neu, IL1R2 + Neu, S100A8/9 + Neu, degranulating CEACAM8 + Neu and mature CD10 + CD16 + Neu frequency with the fraction of proliferative CD4 + T cells, PD-1 expression and CD69 expression. None showed statistically significant correlation (Extended Data Fig. 2f), suggesting that all Neu subsets might be involved in the observed effect of CD66b + Neu on CD4 + T cell proliferation and PD-1 and CD69 expression.
Depletion of arginine by increased arginase-1 activity and upregulation of T cell immune checkpoint PD-1-PD-L1/PD-L2 pathways can modulate the immunosuppressive properties of granulocytic myeloid-derived suppressor cells (G-MDSCs) 12,15 . Inhibition with arginine and an Arg1 inhibitor, or antibodies to PD-L1/PD-L2 applied to the CD66b + Neu-CD4 + T co-cultures, did not reverse the effects on proliferation, PD-1 and CD69 expression compared to co-cultures without inhibitors (Extended Data Fig. 2g). Addition of a prostaglandin EP2 receptor antagonist (TG6-10-1) increased CD69 expression compared to no inhibitor (Extended Data Fig. 2g), whereas the cyclo-oxygenase inhibitor indomethacin or a selective prostaglandin EP4 receptor competitive antagonist (GW-627368) did not (Extended Data Fig. 2g), suggesting effects at the level of pre-formed prostaglandin E2.
Multimodal clustering on 29,336 hematopoietic stem cells (HSCs) retained after excluding progenitor cells identified seven clusters (Fig. 3a). Comparing patients with sepsis to HCs, clusters C3, C4, C5 and C7 were enriched in sepsis, whereas C6 was enriched in HCs ( Fig. 3b and Extended Data Fig. 3j). To understand whether any expanded clusters represented HSCs with a granulopoietic bias, we leveraged Human Cell Atlas (HCA) annotated HSCs and progenitor scRNA-seq bone-marrow data 26 to define genes upregulated in granulopoietic cells, including MPO, RNASE2, ELANE and FKBP2 (Methods and Extended Data Fig. 3k). Clusters C5 and C7 had the highest expression of this gene set (Fig. 3c). Granulopoiesis is driven by transcriptional circuits mediated by CEBP transcription factors, where CEBPA is important during steady-state granulopoiesis (SSG) and CEBPB is a key regulator of emergency granulopoiesis (EG) 27,28 . CEBPA/CEBPB binding motifs were enriched in accessible chromatin sites in both C5 and C7 (Fig. 3d), consistent with upregulated SSG and EG. Comparison of gene expression profiles of C3, C4, C5 or C7 with C6 identified genes, including SETBP1 and FLT3, which are involved in expanded myelopoiesis 29,30 , as upregulated in C5 and C7, respectively (Fig. 3e,f). These data provide evidence for altered granulopoiesis in sepsis that involved specific HSC clusters.  Article https://doi.org/10.1038/s41590-023-01490-5

Altered granulopoiesis drives a sepsis subphenotype
We next investigated whether the altered Neu-granulopoietic profile observed here was related to a previously identified SRS1 subphenotype in patients with sepsis 5,17 . Based on the expression of a seven-gene set (DYRK2, CCNB1IP1, TDRD9, ZAP70, ARL14EP, MDC1 and ADGRE3), the SRS1 subphenotype was assigned to 16 out of 26 patients in the scWB cohort (Extended Data Fig. 4a,b), in agreement with consensus or unsupervised clustering of the pseudobulked scRNA-seq data (Extended Data Fig. 4c-e). IL1R2 + Neu and cycling MK167 + CYP1B1 + Neu were increased in patient samples assigned as SRS1 compared to non-SRS1, whereas mononuclear cells, including CD14 + CD16 − cMo, CD56 + natural killer (NK) and LTB + IL7R + memory CD4 + T cells, were depleted ( Fig. 4a,b). Large numbers of differentially expressed genes were detected in mature CD10 + CD16 + Neu (EBI3, MYOSLID and SLC1A3), S100A8/9 hi Neu (SLC1A3, KLF14 and SGPP2), degranulating CEACAM8 + Neu (SFXN1, HS3ST3B1 and PDE4D) and IL1R2 + Neu (NAIP, IDI1 and   subphenotype. a, Differential cell abundance from whole-blood scRNA-seq UMAP with sampled neighborhoods colored by statistical significance (spatial FDR < 0.05) in patient samples in the scWB cohort (n = 26 patients with sepsis) assigned as SRS1 (n = 16) or non-SRS1 (n = 10). b, Beeswarm plot of differential cell abundance in scWB with cluster labels of neighborhoods depicted and compared for patients assigned as SRS1 and non-SRS1 as in a. c, First two principal components from PCA of pseudobulked Neu states colored by SRS assignment. d, Percentage of PD-1 + (top) and CD69 + (bottom) CD4 + T cells after 72-96 h of co-culture with CD66b + Neu isolated from samples assigned as SRS1 (n = 6) and non-SRS1 (n = 13). Box plots denote minimum and maximum with whiskers and bottom quartile, median and upper quartile with the box. e, Assessment of phagocytosis in CD66b + Neu isolated from samples assigned SRS1 (n = 15) and non-SRS1 (n = 14) and incubated with pHrodo Green E. coli bioparticles stained with 7-AAD, CD66b-AF700 and Siglec-8-APC. Box plots denote minimum and maximum with whiskers and bottom quartile, median and upper quartile with the box. Functional assays tested with two-sided Wilcoxon rank-sum tests. *P < 0.05.

Article
https://doi.org/10.1038/s41590-023-01490-5 CLEC4D) in samples assigned as SRS1 compared to non-SRS1 (Extended Data Fig. 5a), with a corresponding separation by SRS status on principal-component analysis (PCA) for these populations (Fig. 4c). By contrast, there were minimal differences in mononuclear cell subsets between SRS1 and non-SRS1 (Extended Data Fig. 5b), consistent with SRS groupings being driven by Neu. Cell surface expression of IL-1R2, measured by flow cytometry, showed a moderate correlation with expression of IL1R2 in CD66b + Neu isolated from patient samples assigned SRS1 (Extended Data Fig. 5c). In co-cultures, SRS1 CD66b + Neu suppressed CD4 + T cell activation more than non-SRS1 CD66b + Neu (Fig. 4d), but not CD4 + T cell proliferation (Extended Data Fig. 5d). Moreover, SRS1 CD66b + Neu displayed reduced phagocytosis compared to non-SRS1 (Fig. 4e), which was not restored by the addition of granulocyte-macrophage colony-stimulating factor (GM-CSF) (Extended Data Fig. 5e). These data indicated that SRS1 represented an immunosuppressed state driven, at least in part, by Neu dysfunction.
To investigate this further in a multimodal validation (mmV) cohort, we reanalyzed whole-blood bulk RNA-seq and mass cytometry (CyTOF) immunophenotyping data 31 in 42 samples from 36 individuals with all-cause sepsis and 11 age-and sex-matched HCs (Supplementary Table 1). SRS assignment of the 36 patients recapitulated higher early mortality in SRS1 compared to non-SRS1 patients (Extended Data Fig. 6a,b). We identified eight Neu clusters in the CyTOF dataset (Extended Data Fig. 6c-f). More immature CD64 + CD10 lo CD16 lo CD15 lo Neu and CD71 hi CD38 + Ki-67 + pro-Neu or CD71 lo Ki-67 + pre-Neu were detected in patient samples assigned SRS1 compared to non-SRS1 (Extended Data Fig. 6g,h). Pseudotime trajectory analysis, which arranged cells in a progression of sequential maturation stages, demonstrated over-representation of SRS1 samples earlier in the trajectory (Extended Data Fig. 6i-l). PCA of CyTOF data showed separation of Neu subsets, but not mononuclear cell compartments, between samples assigned as SRS1 compared to non-SRS1 (Extended Data Fig. 6m), indicating changes in Neu subsets drove the patient SRS grouping.
To further investigate drivers of the SRS subphenotype, we reduced mmV RNA-seq dimensionality to 33 gene modules using weighted gene coexpression network analysis (WGCNA) and identified 13 differentially expressed modules between SRS1 and non-SRS1 assigned samples (FDR < 0.01) (Extended Data Fig. 7a). The SRS1 upregulated modules were enriched for gene expression signatures of differentiating Neu 32 . Module 10 (bone-marrow Neu and stage 1 differentiating Neu gene sets, including ALOX5, CYBB and LCN2) (Extended Data Fig. 7b) correlated with immature CD64 + CD10 lo CD16 lo CD15 lo Neu frequency (Extended Data Fig. 7c) and a gene expression signature defining IL1R2 + Neu (in the scWB cohort) including IL1R2, PFKFB2 and RETN genes (Extended Data Fig. 7d and Methods). To further validate enrichment of IL1R2 + Neu in SRS1, we performed cell type and cell state deconvolution on an independent total leukocyte microarray dataset 5,17,18 consisting of 542 patients with all-cause sepsis, using the scWB single-cell multiomics dataset as a reference. We observed that IL1R2 + Neu and cycling MK167 + CYP1B1 + Neu were increased in patient samples assigned as SRS1 compared to non-SRS1 (Extended Data Fig. 7e).
Next, we performed multiomics factor analysis (MOFA) to integrate cell cluster (CyTOF), gene module (RNA-seq) and plasma proteomic data (Extended Data Fig. 7f and Methods). Of all factors, 1 and 2 were most divergent between patient samples assigned as SRS1 compared to non-SRS1 (Extended Data Fig. 7g-k). The SRS1 direction in both factors was driven by immature Neu and progenitor Neu cell abundance (including CD71 hi CD38 + Ki-67 + pro-Neu, CD71 lo Ki-67 + pre-Neu and CD64 + CD10 lo CD16 lo CD15 lo immature Neu), together with SRS1 upregulated gene modules and accounted for most variance in the cell and gene module datasets, but not proteomics (Extended Data Fig. 7i). Factor 5 was most strongly related to differences in plasma protein abundance, but showed no difference across SRS (Extended Data Fig. 7g-i), consistent with gene modules and cell clusters, rather than plasma proteins, being associated with SRS. Overall, the results show that differences in immature Neu populations, assayed using different modalities and Neu-granulopoietic dysfunction, were enriched in the SRS1 patient subphenotype.
Given the role of STAT3 in EG and enrichment for STAT3 in SRSq-correlated Neu-GEP, we investigated whether cytokines known to induce granulopoiesis and signal through STAT3 differed in plasma abundance (mmV cohort). Granulocyte colony-stimulating factor (G-CSF) and IL-6 were elevated in patient samples assigned as SRS1 compared to non-SRS1, whereas macrophage colony-stimulating factor (M-CSF) and GM-CSF showed no statistically significant difference (Fig. 5d). Consistent with G-CSF priming of Neu in cancer increasing NETosis 33 , we found that CD66b + Neu from patient samples assigned as SRS1 underwent more NETosis than non-SRS1 samples on live-cell imaging for DNA-bound Cytotox Green reagent when stimulated with phorbol 12-myrisate 13-acetate (PMA) (Fig. 5e,f). These data demonstrated elevated STAT3-driven GEP in different SRS1-correlated Neu subsets and increased levels of circulating cytokines that signal through STAT3 and control granulopoiesis in SRS1 compared to non-SRS1 patients.

STAT3 processes in HSC cluster drive EG in SRS1
To test whether the SRS1 patient subphenotype represented a state of maladaptive granulopoiesis that involved heightened STAT3-mediated EG, we investigated whether there was any difference in circulating HSCs between SRS1 and non-SRS-1 subphenotypes, specifically the cell clusters differentially associated with granulopoiesis (C5 and C7), in 15 patients with sepsis from the single-cell hematopoietic stem and progenitor cell (scHSPC) atlas cohort. Cluster C5, but not C7, was enriched in patient samples assigned SRS1 compared to non-SRS1 (Fig. 6a,b and Extended Data Fig. 8a). We characterized the differentially accessible regions (DARs) of chromatin that defined the clusters and found 184 and 718 DAR in C5 and C7, respectively compared to other clusters (FDR < 0.05, fold change (FC) > 1.5) (Extended Data Fig. 8b,c and Supplementary Table 5). The chromatin profiles for C5 and C7 were enriched for publicly available myeloid progenitor cell chromatin profiles (C5, HSC multipotential progenitor; C7, megakaryocyte progenitor) (Extended Data Fig. 8d,e), suggesting their myelopoietic bias. CEBP motifs were identified within the DARs in both C5 and C7 (Extended Data Fig. 8f,g), whereas STAT motifs were enriched in C5, but not C7 DARs (Extended Data Fig. 8f,g). To differentiate which transcription factors governed the identities of clusters C5 and C7, we overlapped Article https://doi.org/10.1038/s41590-023-01490-5 Cohort 3 (module 10) S100A8/9 hi Neu (program 8)

Article
https://doi.org/10.1038/s41590-023-01490-5 Differential HSC abundance between patients with sepsis assigned as SRS1 (n = 6) or non-SRS1 (n = 9) assessed using density-based distribution UMAP visualization (a) and the corresponding beeswarm plots with sampled neighborhoods colored by statistically significant enrichment (spatial FDR < 0.05) (b). c, chromVAR transcription factor motif z score deviation for STAT3 in HSC from patients with sepsis (n = 15  DARs with public chromatin immunoprecipitation (ChIP)-seq datasets to match transcription factor occupancy profiles. Enrichment of binding profiles for both CEBPA and CEBPB was detected in C7 (Extended Data Fig. 8h), whereas only CEBPB was enriched in C5, indicating that C5 was more biased toward EG while C5, but not C7, exhibited overlap with STAT3 binding profiles (Fig. 6c,d).
To identify transcription factors relevant to HSC clusters, we integrated RNA and chromatin accessibility from the scHSPC data. We constructed supervised pseudotemporal trajectories from cluster C6 (enriched in HCs) to C5 (enriched in SRS1) and C7 (enriched in sepsis, but not SRS1) ( Fig. 6e and Extended Data Fig. 9a). STAT3 and CEBPB expression increased along the C6-C5 trajectory (Fig. 6e), but not the   and pediatric septic shock (n = 106). Box plots denote minimum and maximum with whiskers and bottom quartile, median and upper quartile with the box. Twosided Wilcoxon rank-sum test comparing SRS1 and non-SRS1 groups. b,c UMAP (b) and frequency of immature IL1R2 + Neu in samples assigned as SRS1 or non-SRS1 in adult ARDS patients (n = 9) (c) after reference mapping of scRNA-seq data from whole blood to reference single-cell atlas derived from the scWB showing. Box plots denote minimum and maximum with whiskers and bottom quartile, median and upper quartile with the box. d, Heat maps of expression of STAT3 GEPs (CD10 + CD16 + Neu_program_3, S100A8/9 hi Neu_program_8 and IL1R2 + Neu_program_8) in the corresponding Neu subsets from patients with ARDS (n = 9) with violin plots on each side. Two-sided Wilcoxon rank-sum test comparing SRS1 and non-SRS1 groups. e, Violin plots of expression of a combined set of GEPs involving STAT3 (CD10 + CD16 + Neu_program_3, S100A8/9 hi Neu_ program_8 and IL1R2 + Neu_program_8) from scWB Neu subsets in pseudobulked Neu from patients with ARDS (n = 9) assigned as SRS1 or non-SRS1. Two-sided Wilcoxon rank-sum test comparing SRS1 and non-SRS1 groups. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.
To validate the importance of CEBPB-STAT3 and CEBPA in governing the identities of C5 and C7, respectively, transcription factor-mediated gene regulatory networks were constructed with CellOracle from RNA-seq and ATAC-seq HSC data for all HSC clusters. In silico knockout of CEBPA led to C7 loss of identity and transition toward C6, whereas C5 cells were not affected (Extended Data Fig. 9e), while in silico CEBPB knockout disrupted C5 cell differentiation (Extended Data Fig. 9e). In silico overexpression of CEBPA or CEBPB reversed the directions of cell fate transitions in C5 and C7, respectively (Extended Data Fig. 9e), suggesting that C7 was a SSG cluster and C5 was an EG cluster. In silico knockout or overexpression of STAT3 had similar effects as CEBPB (Fig. 6g), suggesting that STAT3 was driving EG. An independent methodology, using scRNA splicing information and vector field analysis of HSCs C5, C6 and C7, showed the same effects on EG and SSG reversal following in silico knockout of CEBPB, CEBPA and STAT3 (Extended Data Fig. 9f-j). These observations established C7 as a CEBPA SSG HSC cluster and C5 as a CEBPB-STAT3 EG cluster, with SRS1 enrichment of C5 highlighting EG in SRS1.
We also analyzed published whole-blood transcriptomic datasets of patients with ARDS 14 (n = 9) or COVID-19 (ref. 13) (n = 8) at a single-cell resolution (scWB-CID) using scWB as a reference for cell type and cell state annotation (Fig. 7b). IL1R2 + Neu were expanded in patients with SRS1 subphenotype compared to non-SRS1 ( Fig. 7c and Extended Data Fig. 10c). The STAT3 GEP (CD10 + CD16 + Neu_program_3, S100A8/9 hi Neu_program_8 and IL1R2 + Neu_program_8) showed enrichment in patient samples assigned as SRS1 compared to non-SRS1 in ARDS (Fig. 7d) and COVID-19 (Extended Data Fig. 10d). Total Neu, as defined by the original authors 13,14 , showed higher expression of the STAT3_combined_program genes in patient samples assigned as SRS1 compared to non-SRS1 in ARDS (Fig. 7e) and COVID-19 (Extended Data Fig. 10e) providing cross-validation. These analyses indicated that the biological basis of the SRS patient subphenotype was independent of infectious organism, source of infection, clinical syndrome and age.

Discussion
Here, we showed that Neu-granulopoietic disturbances in sepsis involved expansion of specific populations of immature Neu, suppression of CD4 + T cells in co-culture and altered granulopoiesis, and demonstrated these features were enriched in a subset of patients (SRS1). These results defined SRS1 as a specific immunocompromised disease endotype.
Our fresh whole-blood single-cell multiomic atlas, with no cellular enrichments or depletions, ensured faithful recapitulation of the sepsis cellular landscape. Previous single-cell -omic profiling focused on PBMCs, identifying an immature, bone-marrow-derived monocyte state (MS1), expanded in and predictive of sepsis 8 . The immature Neu populations defined here, in particular IL1R2 + Neu, exhibit similar gene expression profiles to MS1 cells 8 . Direct comparison between MS1 and IL1R2 + Neu could potentially reveal similar myelopoietic processes leading to their generation. The extent to which mobilization of IL1R2 + Neu may occur elsewhere in systemic inflammation, such as reported in mice 37 and whether IL1R2 + Neu and other Neu subsets may drive immune suppression in these contexts, remains unclear.
While altered myelopoiesis in sepsis has been described in mice 38,39 and modeled in vitro 40 here we presented evidence for amplified granulopoiesis in humans during sepsis and specifically, a dysregulated form of EG in the SRS1 endotype. Our data, together with the reported increased risk of infections in patients with clonal hematopoiesis, including sepsis 41 , triangulate on the bone marrow as foundational for the maladaptive response to infection, with the caveat that our samples derive from circulating HSPCs rather than bone-marrow tissue.
Going forward, it will be important to understand determinants of differential bone-marrow responses to infection and the myelopoietic legacy of severe infection, for example, through trained immunity 42 or hematopoietic exhaustion 43 . Published reports suggest that previous exposures and inflammatory comorbidity may be important in influencing subsequent myelopoietic responses to infection 44 . Meanwhile, our observations of persistent granulocytic alterations in convalescence, both phenotypic and functional, add to the evidence that infectious and inflammatory stimuli have long-lasting myelopoietic and therefore innate immune ramifications 42,45 . CEBPB has a key role in induction of trained immunity in HSCs 46 and STAT3 drives specific immunosuppressive properties in MDSCs 21 . We found that STAT3-CEBPB-driven EG was pathognomonic of SRS1, with STAT3 underpinning the granulopoietic-granulocytic axis, raising the hypothesis that SRS1 represents a state of maladaptive innate immune reprogramming and memory, with a potential opportunity for manipulating STAT3 activation to alleviate sepsisand SRS1-associated immunosuppression. For example, G-CSF and IL-6 show increased expression in SRS1, canonically signal through STAT3 and are key contributors to EG 47,48 with inhibition specifically in patients with SRS1 or high SRSq subphenotypes a possible immunotherapeutic strategy 1 . This is further supported by Mendelian randomization work, where lower IL-6R expression associated with reduced mortality in sepsis 49 and the therapeutic benefit of targeting IL-6 in severe COVID-19 (ref. 50).
Limitations of our study include the extent that our single-cellanalysis patient cohorts are fully representative of the breadth of the sepsis syndrome. Further work is needed to understand the differential immunosuppressive properties and role of prostaglandins in Neu, altered granulopoiesis through study of patient bone marrow and experimental gene manipulation of transcription factors to verify in silico knockouts. Relevant animal model and experimental medicine studies to manipulate candidate therapeutic targets are needed to better understand cytokine inhibition strategies.
Collectively, our work identified a common innate immune and hematopoietic axis that contributes to the maladaptive immune response to infection during sepsis and specifically a poor outcome, Article https://doi.org/10.1038/s41590-023-01490-5 immunocompromised, extreme response SRS1 patient endotype, advancing opportunities for personalized medicine.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41590-023-01490-5.

Fresh whole-blood sample processing for single-cell sequencing
Blood samples were drawn into EDTA tubes (BD Biosciences) and processed within 1 h of collection. Then, 1 ml whole blood was lysed with 9 ml 1× eBioscience red blood cell lysis buffer (Thermo Fisher) twice and 200,000 cells were transferred for antibody staining.
Whole-blood single-cell RNA and cell surface protein profiling scRNA and cell surface protein sequencing was performed with the BD Rhapsody platform (633731/633733, whole transcriptome assay (633801)) using 30 AbSeq antibodies (1 μl per antibody) (Supplementary Table 2). Cells were stained following the manufacturer's recommendations before single-cell capture targeting 6,000 cells per sample. Reverse transcription, complementary DNA amplification and library construction (633801) were performed following the manufacturer's recommendations in six batches. Libraries were sequenced on a NovaSeq6000 (Illumina).

Whole-blood single-cell multiomic analysis
Analysis was performed with v.4.0.0 R and Python 3.8.6.
Preprocessing and quality control. Gene expression data were aligned using STARsolo 53  Normalization, dimensionality reduction and clustering. Data were log normalized and 4,000 highly variable genes (HVGs) were identified using the Seurat vst algorithm (scanpy v.1.7.2).
Multimodal dimensionality reduction was performed with TotalVI (scvi-tools v.0.10.0) 55 on default settings with all 30 proteins and HVGs with each individual sample set as a batch (as all samples were processed separately).
Unsupervised clustering was performed on the 20 TotalVI latent dimensions (Seurat FindNeighbors (k = 30) and FindClusters (smart local moving algorithm)). Clustering resolution was evaluated by cluster neighborhood purity, cluster average silhouette width and a 30-iteration bootstrap to determine cluster stability with respect to sampling noise (bluster v.1.0). Additionally, we inspected top markers per cluster to match with known biology and understand potential value in merging versus splitting clusters. Combining these elements, the default clustering resolution 0.8 was chosen, with one immature neutrophil cluster split into two (MPO+ immature neutrophils/progenitors and PADI4+ immature neutrophils) based on resolution 1.1.

Cell annotation and RNA velocity/trajectory analysis.
Cell clusters were merged for protein marker-based annotation of major known immune cell types at a broad level and kept at the clustering resolution of choice for fine annotation. Lineage assignment was confirmed with SingleR (v.1.6.1) assignments. Fine annotation was conducted by inspecting biologically meaningful gene markers for T cell and neutrophil populations (Supplementary Table 3).
RNA velocity and partition-based graph abstraction analysis of the neutrophils (without degranulating or apoptosing neutrophils) was performed with scVelo (v.0.2.3, default settings of stochastic model) 56 . Moments were estimated using TotalVI reduced dimensions.
scRNA-seq pseudobulk unsupervised analysis (consensus clustering, PCA, hierarchical clustering) and SRS assignment/scoring. Gene expression was aggregated per individual for all cells and normalized (EdgeR trimmed mean of M-values; TMM). The top 10% most variable genes by mean absolute deviation were taken for unsupervised analysis.
The relevant genes (DYRK2, CCNB1IP1, TDRD9, ZAP70, ARL14EP, MDC1 and ADGRE3) were used for sepsis sample SRS assignment and SRSq score calculation by SepstratifieR (v.0.0.0.9) 57 . Differential abundance analysis. Cell type/state DA across conditions was identified by sampling neighborhoods of cells from a k-nearest neighbors (k-NN) graph and looking for enrichment of either condition in each neighborhood as implemented in MiloR 58 . The 20 batch-corrected latent dimensions from TotalVI were used for MiloR (v.0.99.19) k-NN graph construction (k = 30) and neighborhood indexing (proportion = 0.1). DA testing was performed with generalized linear models, including age and sex as covariates (neighborhoods significant if spatial corrected FDR < 0.05).

DGE analysis.
For pseudobulk DGE, gene expression was aggregated per individual and per cell type/state into pseudobulks. Genes were filtered per pseudobulk based on minimum expression of n counts in at least X samples, where X was the smallest comparator group and n was defined for each pseudobulk based on histograms of logged count distributions. Pseudobulks were normalized by the EdgeR TMM method (v.3.30.3) 59 . DGE was performed with generalized linear models as implemented in EdgeR with age, sex and sequencing batch included in the model.
Consensus DGE for acute sepsis versus HC IL1R2+ immature neutrophils (as cells were too sparse to run pseudobulk DGE) and convalescent versus HC total neutrophils was performed as described in ref. 60. Genes were filtered for those that showed reproducible change in the same direction in a minimum of six samples (smallest comparator group size, HC = 6 samples).
Consensus non-negative matrix factorization of neutrophil gene expression. The 1,000 HVGs for each neutrophil state were selected (Seurat SelectIntegrationFeatures) followed by cNMF (v.1.2) as previously described (100 factorizations) 61 for 4-15 GEPs. The final GEP number per neutrophil state was chosen based on a tradeoff between stability versus error as recommended by the original authors. Mean GEP usage per sample for each neutrophil state was correlated with sample SRSq scores with FDR adjustment. The top 50 genes per GEP were taken for downstream analysis.
Transcription factor prediction analysis. GEPs positively correlated with SRSq were used for TF prediction analysis with the Cytoscape (v.3.9.1) plugin iRegulon (v.1.3) 62 (default settings). The same analysis was performed for the top 1% genes correlating with module 10 eigengene from mmV cohort RNA-seq WGCNA (below). IL1R2 + neutrophil defining gene set. Genes specific to IL1R2+ neutrophils were identified with Wilcoxon tests (Seurat FindMarkers) for sepsis IL1R2 + neutrophils, filtering to only retain genes with adjusted P value <0.01 and average log 2 FC > 1.

Neutrophil functional assays
Isolation of neutrophils from whole blood. Neutrophils were isolated from 5-10 ml whole blood from EDTA Vacutainer tubes (BD Biosciences) using EasySep HLA Chimerism Whole Blood CD66b positive selection kit (StemCell) following the manufacturer's instructions.
Phagocytosis assay. Neutrophils were incubated for 20 min in complete medium with or without (fluorescence minus one) pHrodo Green E. Coli bioparticles (Invitrogen) (1 neutrophil:10 bioparticles). Cells were washed and stained for surface markers for 30 min with 7-AAD, CD66b-AF700 (G10F5) and Siglec-8-APC (7C9) antibodies from BioLegend and acquired with a BD LSRFortessa X-20 analyzer. The phagocytosis median fluorescence intensity for single 7-AAD − CD66b + Siglec-8 − cells was determined by subtracting the median fluorescence intensity of the fluorescence minus one control sample.
After 72-96 h of co-culture, cells were stained for surface markers followed by annexin V and 7-AAD to exclude dead/apoptosing cells. Samples were analyzed using the BD LSRFortessa X-20. All antibodies were purchased from BioLegend unless otherwise stated: CD3-APC (UCHT1), CD4-BUV395 (BD Biosciences, SK3), CD66b-AF700 (G10F5), PD-1-PE (NAT105), CD69-PECy7 (FN50), 7-AAD, annexin V − FITC. T cells were gated as singlet CD66b − CD4 + CD3 + annexin V − 7-AAD − cells. Proliferation analysis by dye dilution was established using the non-bead-stimulated T cells cultured without neutrophils as the baseline proliferative fraction of cells in each sample. The proliferative fraction and percentage of cells expressing PD-1 and CD69 were calculated as a percentage relative to the anti-CD3/28 bead-stimulated, no co-culture T cells control to account for donor variation.

NETosis assay.
Neutrophils were plated at 20,000 cells per well (100 μl per well) in Ham's F-12K medium (Gibco) using 96-well flat-bottom plates coated with 0.01% poly-l-ornithine solution (Sigma). Then, 250 nM IncuCyte Cytotox Green Dye was added to measure NETosis. Cells were treated with 100 nM PMA as per the manufacturer's recommendations then imaged with the IncuCyte Live-Cell Analysis System for up to 4 h.

PBMC isolation and cryopreservation
PBMCs were isolated from sepsis (acute and convalescent) and HC whole blood using density gradient centrifugation with Leucosep tubes (Greiner) and lymphoprep (StemCell) and cryopreserved in 10% dimethylsulfoxide (Cell Signaling Technology).

CD34 + hematopoietic stem and progenitor cell enrichment
PBMCs from 15 patients with acute sepsis, 7 age and sex-matched HCs and 8 convalescent sepsis samples were used for CD34 + HSPC isolation.

Nature Immunology
Article https://doi.org/10.1038/s41590-023-01490-5 Three convalescent sepsis samples were removed after preprocessing and demultiplexing (below) (samples were judged as too long since the acute episode (>6 months past hospital discharge) by retrospective clinical evaluation). The 30 samples were processed in six batches, with each batch containing at least one sample from each comparator group.

HSPC single-cell multiomics
scRNA and scATAC-seq on CD34 + HSPCs was performed by isolating nuclei of flow-sorted HSPCs. An equal number of cells per sample was pooled for each batch. Cell lysis and nuclei extraction were conducted following the low input workflow within the 10x Genomics Demonstrated Protocol (CG000365 Rev B). Nuclei were transposed, captured and RNA and ATAC library preparation (1000285, 10x Genomics) was performed as per the 10x Genomics manufacturer's protocol. RNA libraries were sequenced on a NovaSeq6000 (Illumina) and ATAC libraries were sequenced on a NextSeq500 (Illumina).
Quality control. Cells expressing <100 or >6,000 genes, >25,000 UMIs or with a log 10 (UMI per gene) <0.8 were removed. Cells with TSS enrichment <7 and <1,000 unique fragments were also filtered out. Genes expressed in <10 cells or <3 in total count were removed. This left 46,782 cells (median TSS enrichment of 14.6, median fragment count of 13,821) and 26,660 genes.
Cell type/state annotations and non-HSC filtering. We assigned identity of HSPCs and contaminating cells by mapping our scRNA-seq data to two healthy donor reference bone-marrow mononuclear cell scRNA-seq datasets from ref. 24 (reference dataset 1, supervised PCA) and ref. 25 (reference dataset 2, PCA) with Seurat (v.4.0). Cells with non-HSPC labels from either mapping were filtered out, leaving 46,156 HSPCs.
Multimodal dimensionality reduction for HSPCs was performed (iterativeLSI) for scRNA (gene expression matrix) and scATAC (tile matrix). Both modalities were batch-corrected (Harmony) and combined (addCombinedDims) before clustering on the combined dimensions (resolution of 0.4, smart local moving algorithm). HSPC identity was finalized by assigning the majority RNA mapping identity from reference dataset 1 per cluster. Progenitor cells were thus filtered out and only 29,336 HSCs retained for downstream analysis.
Progenitor cell differential abundance analysis. Hematopoietic progenitor cells from the above filtering step (n = 16,820) were compared across sepsis and HC conditions by two-sided Wilcoxon rank-sum tests with FDR adjustment for multiple testing to look for differential abundance.
HSC peak calling and reclustering. Each individual sample was pseudobulked for peak calling with MACS2 and iterative peak overlapping removal 65 within ArchR (minimum cells = number of cells in the sample with the fewest cells, minimum replicates = 5 (lowest sample group size that is convalescent sepsis samples)). We re-performed dimensionality reduction and batch correction on the gene expression and peak matrices before reclustering HSCs on combined dimensions (smart local moving algorithm, resolution = 0.2 to avoid overclustering and challenging interpretability given the purity of the cell type already). Cluster C2 was derived from only one single sample.
HSC differential abundance analysis. HSC cluster DA changes across conditions were conducted in MiloR (v.1.2.0). The 60 batch-corrected reduced dimensions (30 each of scRNA/scATAC) were used in MiloR k-NN graph construction (k = 20) and neighborhood indexing (proportion = 0.3). DA testing was performed with generalized linear models, including age and sex as covariates (neighborhoods were significant if spatially corrected FDR < 0.05).

Granulopoiesis gene set enrichment scoring for single cells.
We reanalyzed HCA HSPC scRNA-seq data, utilizing annotations from Hay et al. 26 to contrast all granulocytic progenitors with HSCs, multilineage progenitors, monocyte/DC progenitors, lymphoid-primed multipotent progenitors and megakaryocyte/erythroid progenitors to define a granulopoietic gene set. We scored our HSCs with this gene set (Seurat AddModuleScore) and compared differential gene set activity across HSC clusters with a Kruskal-Wallis test and post hoc Dunn's test with FDR adjustment.

DGE and chromatin accessibility analysis.
Genes defining the C5 and C7 HSC clusters (versus C6) were identified with Wilcoxon tests (get-MarkerFeatures) with thresholds of FDR < 0.05 and FC > 1.5. Genomic regions defining the C5 and C7 HSC clusters were identified with Wilcoxon tests (getMarkerFeatures) with thresholds of FDR < 0.05 and FC > 1.5.

Bulk ATAC/TF motif enrichment analysis and ChIP-seq dataset overlapping.
Differentially open peak regions for clusters C5 and C7 were taken for bulk ATAC/TF motif enrichment analyses (peakAnno Enrichment). We also input the differentially open peak regions into the cistromeDB toolkit (http://dbtoolkit.cistrome.org/) 66 (top 1,000 peaks according to peak enrichment used).
Gene regulatory network construction and in silico knockout/ overexpression analyses. We built gene regulatory networks (GRNs) with CellOracle (v.10.10) 67 following tutorials from https://morris-lab. github.io/CellOracle.documentation/. In brief, HSPC scATAC data were used to build the base GRN via Cicero 68 , following which, scRNA data were used to prune the GRN though a Bagging Ridge model. Instead of recomputing principal components for the scRNA data for GRN construction, we instead used the batch-corrected joint ATAC and RNA LSI dimensions. Effects on gene expression of HSC clusters C5/6/7 when perturbing specific TFs (CEBPB, CEBPA and STAT3) were simulated by setting expression of the TF either to 0 or double the value of its maximum expression, with CellOracle then predicting based on the new simulated gene expression changes the trajectory of cellular transition. Article https://doi.org/10.1038/s41590-023-01490-5 RNA velocity and vector field analyses. We recapitulated expression dynamics vector fields with dynamo (v.1.1.0 (ref. 69)) following tutorials from https://dynamo-release.readthedocs.io/en/latest/. In brief, HSPC scRNA raw sequencing data that are FASTQ files were reprocessed by STARsolo (v.2.7.9a) (GRCh38) and spliced and unspliced counts produced for the already identified HSPCs. RNA velocity, acceleration and curvature were then calculated for clusters C5/6/7. Effects on RNA velocity on HSC clusters C5/6/7 when perturbing specific TFs (CEBPB, CEBPA and STAT3) were simulated by setting expression of the TF to 0. The effect of a change in expression of these TFs on the regulation (activation versus inhibition) of each other was also analyzed by calculating the RNA Jacobian for the HSCs.

SRS assignment for neutrophil functional assay samples and scHSPC samples not part of scWB or mmV cohorts
Whole blood was sampled into Tempus tubes (Thermo Fisher) and stored at −80 °C until RNA extraction with Norgen Preserved Blood RNA Purification kit I (43400) according to the manufacturer's instructions. RNA was used for quantitative PCR with reverse transcription for the seven SRS genes and Cq values input into SepstratifieR for SRS assignment as previously described 57 .

Bulk RNA-seq analysis
Analysis was performed with v.4.0.0 R. The EdgeR TMM normalized count matrix was obtained from the COMBAT consortium 31 and SRS assignment was conducted as for scWB samples.
Pathway analysis. Pathway analysis was performed against neutrophil gene sets from the HCA 26 and signatures of differentiating neutrophils from mouse single-cell transcriptomics 32 (ClusterProfiler enricher). Gene set enrichment analysis for scWB IL1R2 + neutrophil defining genes was performed via the gene set enrichment analysis function against all genes in module 10 ranked by their correlation with the module eigengene.
Preprocessing. Normalized, debarcoded, bead and doublet-cleaned FCS files were obtained from the COMBAT consortium 31 . To avoid biases due to highly varying cell numbers per sample, a maximum of 50,000 cells per sample was taken. Harmony (v.1.0) 70 was used for batch correction 71 .
Clustering, trajectory inference and differential abundance analyses. Clustering and consensus clustering were performed with the CATALYST package (v.1.14.0) 72,73 with all cells and 41 markers. The 50 metaclusters were merged following manual annotation. As the 41-marker panel was myeloid focused, we annotated five monocyte clusters, before subsetting the neutrophils and reclustering with 17 selected markers. Thirty neutrophil metaclusters were manually annotated and merged based on median marker expression 74 into eight subsets, leaving a final total of 22 clusters at the finest level of annotation (14 non-neutrophil and eight neutrophil).
For trajectory analysis, neutrophils were downsampled to 15,000 cells per comparator group for CATALYST diffusion map dimensionality reduction followed by principal curve fitting (Slingshot v.1.8.0, default settings) 75 .
PCA of cluster proportions was performed after quantile normalization. Cluster DA analysis with FDR adjustment was performed with generalized linear mixed models (glmmPQL function, MASS package, v.7.3-51.6) to account for repeated sampling in six patients, with batch included as an additional fixed effect and patient a random effect. An extra random effect term per sample was included to model overdispersion in proportions seen in high dimensional cytometry data 73 .

Multimodal analysis by MOFA+
We integrated mmV bulk RNA-seq and CyTOF data with MOFA+ (v.1.0.1) 76 , including 105 plasma protein measurements (timsTOF mass spectrometry) from the same patient samples 31 as a computational negative control, reasoning that the latent spaces linking cell clusters to gene expression profiles should not show an overly strong contribution from the plasma proteome, which reflects contributions from various tissues, whereas the RNA was derived solely from leukocytes.

Cytokine analysis
Concentrations of plasma analytes from the Luminex assay were obtained from the COMBAT consortium 31 .

brWB-CID bulk transcriptomic analysis
A signature matrix for CIBERSORTx 77 was constructed using our scWB single-cell dataset, with each finely annotated population except apoptosing cells downsampled to 100 cells per population. The signature matrix was then created via the Create Signature Matrix analysis module with min.expression = 0.25, replicates = 100 and sampling = 0.5.
SRS was assigned for the bulk transcriptomic data of sample set 3 as described for cohorts 1 and 2. Cell fractions of the 1,578 whole-blood bulk transcriptomic samples of sample set 3 were then estimated with the Impute Cell Fractions analysis module using the single-cell reference matrix, with batch correction S-mode enabled quantile normalization disabled for RNA-seq datasets. DGE was performed across SRS for each cohort and SRS1 upregulated genes were tested for enrichment of the union of genes of all three STAT3 GEPs from scWB neutrophil cNMF.

Statistics and reproducibility
In addition to the above details of statistical methods, we ensured analytical rigor with the following procedures.
For scWB, we estimated from existing data 5,18 that the smallest SRS group would be one-third of the recruited patients with sepsis; we opted to recruit 26 patients for a minimum of 8 patients in the smallest group. No sample size calculations were performed. Cells from all study participants were used to determine cell states. For scHSPC, three convalescent sepsis samples were excluded (above), after which all 27 samples were used for all analysis. Sample size was determined to be adequate based on the degree and consistency of differences between groups.
For scWB and scHSPC, sequencing batches were prepared such that each batch contained samples from all comparator groups (acute sepsis, HCs, after CS and convalescent sepsis). SRS groups were not known at the time of patient recruitment and data generation. Cells from all samples in scWB, mmV CyTOF and scHSPC were analyzed blind to which patient they originated from to define the varying cell states.
Data distributions were assumed to be normal for (generalized) linear models but this was not formally tested. Box plots denote minimum and maximum with whiskers and bottom quartile, median and upper quartile with the box.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
Raw data for whole-blood single-cell sequencing and HSC single-cell sequencing datasets are deposited on the European Genome-phenome Archive (EGAS00001006283) and derived data are at Zenodo (https:// doi.org/10.5281/zenodo.7723202). For sequence-level raw datasets deposited at the European Genome-phenome Archive, access is managed by a Data Access Committee.

Code availability
Code used for every algorithm followed in data processing and analysis is fully referenced within the specific Methods sections.