Our study design was previously reported by Darden et al. (14). To summarize, this prospective, observational cohort study was registered with clinicaltrials.gov (NCT02276417) and conducted at a tertiary care, academic research hospital. The objective of the study was to better understand the myeloid response (specifically blood MDSCs) to acute sepsis, and to identify transcriptomic differences in sepsis patients who rapidly recover versus those who develop CCI. Our hypothesis was that differences in the myeloid transcriptomic landscape could explain why some sepsis patients rapidly recover while others develop CCI.

Sepsis designation occurred through an electronic medical record-based Modified Early Warning Signs-Sepsis Recognition System (MEWS-SRS), which uses white blood cell count, heart rate, respiration rate, blood pressure, and mental status to identify patients at risk for sepsis. All patients with sepsis were treated with early goal-directed fluid administration, initiation of broad-spectrum antibiotics, and vasopressor administration if appropriate. Antibiotic treatment was tailored to culture results and antibiotic resistance information.

Inclusion criteria:

Exclusion criteria:

CCI was defined as ICU length of stay >13 days with persistent organ dysfunction as measured by the Sequential Organ Failure Assessment (SOFA) score. Patients were also designated CCI with <14 days ICU length of stay if they were transferred to another hospital, or discharged to a long-term acute care facility or hospice with evidence of persistent organ dysfunction (4, 16).

Whole blood samples were collected at the following time points: 4 ± 1 days and 14-21 days after sepsis (16). For the former, we enrolled 4 patients diagnosed with septic shock (1) in order to guarantee a strong host response and transcriptomic alterations in circulating leukocytes. Interestingly, only half of these septic shock patients went on to develop CCI (17, 18). Samples from 5 patients who developed CCI and 4 patients who rapidly recovered after sepsis were obtained from additional patients between 14-21 days after sepsis diagnosis. We previously determined this time point to be key to MDSC differentiation as well as distinguishing CCI from rapid recovery (4, 16). In addition, whole blood was collected from 12 healthy subjects (1, 2). The proportion of men and women did not differ between sepsis patients and healthy subjects ( Table 1 ). Healthy subjects trended towards being younger than sepsis patients, although they still met the criteria of middle age (>45 years), encompassing patients who have poor outcomes after sepsis compared to younger cohorts (17, 18). Healthy subjects and septic patients had similar underlying comorbidities (most commonly hypertension, chronic obstructive pulmonary disease, and diabetes mellitus).

Each blood collection underwent two enrichment procedures. Peripheral blood mononuclear cells (PBMCs) were collected from half of each sample using Ficoll-Paque™ PLUS (Millipore Sigma, St. Louis, MO) and density gradient centrifugation. Myeloid cells were collected using the RosetteSep™ HLA Myeloid Cell Enrichment Kit (STEMCELL Technologies, Vancouver). A 1:3 mixture of enriched PBMCs: myeloid cells was created in order to adequately analyze the small target population (MDSCs, especially in healthy subjects) while also allowing for characterization of other important immune cell populations using CITE-seq.

We used CITE-seq to analyze single-cell transcriptomic profiles of MDSCs in blood from healthy subjects (n=12) and surgical sepsis patients at 4 ± 1 (n=4) and 14-21 (n=9) days after sepsis onset (65). Septic patients at 14-21 days were further divided based on their clinical outcomes at time of sampling, defined as either ‘rapid recovery’ (n=4) or development of CCI (n=5). Sex, age, and BMI were similar between cohorts ( Table 1 ).

Similar to flow cytometry phenotyping, CITE-seq employs conventional cell surface markers for myeloid cell subpopulations to define E-MDSCs (Lin⁻HLADR^low/− CD33⁺CD11b⁺CD14⁻CD15⁻CD66b⁻), PMN-MDSCs (Lin⁻CD33⁺CD11b⁺CD14⁻ and CD15⁺ or CD66b⁺), and M-MDSCs (Lin⁻HLADR^low/−CD33⁺CD11b⁺CD14⁺CD15⁻CD66b⁻), as well as CD14⁺CD16⁻ (classical) and CD14^dimCD16⁺ (non-classical) monocytes (while removing platelets, erythrocytes, HSPCs, γδ T cells, and innate lymphoid cells). This is consistent with the classical monolithic view of myeloid differentiation described by Hegde ( Figure 1A ) (13). Historically, flow cytometry classification of MDSCs is performed directly on isolated PBMCs (4, 16) and our analysis revealed that PMN-MDSCs made up the majority of MDSCs in isolated PBMCs of representative septic patients, consistent with prior literature ( Table 2 ) (4, 16).

We then evaluated cell proportions using transcriptomics with confirmation via flow cytometry. Myeloid cell enrichment was necessary in this step in order to detect MDSCs in healthy subject samples during CITE-seq as they are a relatively rare population (see Methods Section entitled “Human Blood Collection and Sample Preparation”). Both single-cell transcriptomics and flow cytometry revealed an overall increase in total MDSCs acutely after sepsis ( Figure 1B and Table 3 ). We plotted the cells via UMAP based on timepoint after sepsis ( Figure 1C ) and myeloid cell subtype ( Figure 1D ), which revealed heterogeneity of these cells when analyzing their single cell transcriptomes. The classification of MDSCs based on cell surface phenotypes is somewhat dependent on the method of analysis, and as the myeloid enrichment kit (STEMCELL) uses CD33 (and CD33 is expressed on all MDSCs), this was our first inclination that an alternative method of classifying cells by subpopulation would be necessary.

Next, we performed pseudobulk differential gene expression between septic patients at day 4 and days 14-21 post-sepsis diagnosis compared to healthy subjects to assess possible differences among septic groups. Dramatic differences in gene expression within MDSC subpopulations (specifically PMN- and M-MDSCs, which were most abundant) were observed that varied over time ( Figure 2 ). Considering PMN-MDSCs, by comparing each differential expression test performed against healthy subjects, we found more extreme fold-changes in late sepsis patients who developed CCI compared to acutely septic patients ( Figure 2A , left panel). Conversely, gene expression for late sepsis patients who rapidly recovered returned towards that seen in healthy subjects when compared to both acute sepsis and late sepsis with CCI. Gene expression for rapid recovery and CCI patients compared to healthy subjects also tended to diverge ( Figure 2A , middle and right panels). Overall, 52 genes were differentially expressed in PMN-MDSCs from acutely septic patients; however, only three of these genes were also significantly differentially expressed in patients who rapidly recovered and those who developed CCI ( Figure 2B ; Supplementary File 1 ). The ontology of transcriptional differences among septic patients at different time points also illustrated the heterogeneity of the PMN-MDSC response over time ( Figure 2C ).

The greatest differences in the magnitude of M-MDSC gene expression from healthy subjects compared to septic patients occurred during acute sepsis ( Figure 2D ). In this cohort, 601 genes were differentially expressed in M-MDSCs, and only 31 of these were also significantly differentially expressed in late sepsis patients who rapidly recovered and those who developed CCI ( Figure 2E ). Gene ontology analysis among M-MDSCs revealed that patients experiencing rapid recovery had over-expression of kinase agents versus oxoacid metabolism when compared to sepsis patients with CCI ( Figure 2F ).

In summary, transcriptomic analysis comparing healthy subjects, patients day 4 post-sepsis, and patients at days 14-21 post-sepsis (combining those who rapidly recovered with those who developed CCI) revealed significant heterogeneity in MDSC transcriptomics. In addition, lymphocyte suppressive activity, specifically suppression of T cell proliferation and T cell cytokine/chemokine production, varied between cohorts (see below). MDSCs from both time points after sepsis were dissimilar when comparing their gene expression profiles and significantly enriched biological processes from gene ontology.

Our laboratory has previously used cell sorting and subsequent T-cell suppression assays to demonstrate the immunosuppressive capacity of total MDSCs from septic patients (4, 16). Thus, we set out to verify that functionally active PMN- and M-MDSCs were indeed present in the isolated PBMCs of septic individuals, as defined by Gabrilovich et al. (11). Surprisingly, we discovered unexpected cell types in our flow cytometry samples that were supported by single-cell analysis. Historically in the cancer literature, CD15 positivity is used to distinguish PMN-MDSCs from M-MDSCs (13). However, in representative septic patients, after identifying CD11b⁺CD33⁺ cells ( Figure S1A ) and isolating HLA-DR^-/low cells to obtain the total MDSC population ( Figure S1B ), CD15 was unable to clearly separate MDSC subpopulations ( Figure S1C , bottom panel). Alternatively, CD66b (CEACAM8; a granulocytic marker) was able to better delineate MDSC subpopulations ( Figure S1C , top panel) and, thus, was selected to distinguish PMN-MDSCs from PBMCs for functional analysis ( Figure S1D ) (10, 11, 19, 20, 66). We undertook bulk CD66b⁺ cell isolation (STEMCELL Technologies, Vancouver); however, although CD14^-CD15⁺ PMN-MDSCs enrichment was achieved, further analysis of the CD66b⁺-isolated cells demonstrated distinct CD66b^low and CD66b^high populations ( Figure S2A ). Thus, we had enriched CD66b^lowCD14⁺CD15^- M-MDSCs in our gating strategy which was meant to only contain PMN-MDSCs (CD66b^high) ( Figure S2B ).

Functionally, the CD66b⁺ cells isolated from septic patient PBMCs suppressed either CD4⁺/CD8⁺ T-lymphocyte cytokine/chemokine production ( Figure S3 ) or lymphocyte proliferation of host CD8⁺ T-lymphocytes ( Figure S4 ), thereby meeting the criteria of MDSCs. CD66b⁺ MDSCs from septic but not healthy subjects altered T-cell cytokine production in response to anti-CD3/CD28 treatment, including IFN-γ, IL-2, IL-4, IL-10, and IL-17 ( Figure S3 ). Cytokines which were analyzed which did not exhibit CD66b⁺ inhibition in acutely septic patients include IL-12, IL-23, and TGF-β. CD66b⁺ cells isolated from sepsis patients 14-21 days after infection were capable of significantly suppressing CD8⁺ T-lymphocyte proliferation in response to CD3/CD28 stimulation ( Figure S4 ). Although we did not see suppression of CD4+ T lymphocyte proliferation by CD66b+ cells, we did see a significant decrease in stimulated T-lymphocyte production of IFN-gamma at days 14-21 compared to day 4 (Figure S5). This indicates that T lymphocytes are incapable of maintaining IFN-gamma production 2-3 weeks after sepsis (similar to what has been previously reported) (67), and that MDSCs at this time point may not be able to further suppress this aspect of CD4+ T lymphocyte function (22, 67).

After our initial steps demonstrated that identification of MDSC subsets based on cell surface markers was potentially problematic in sepsis, we transitioned to cell classification via gene expression for the remainder of our analysis. All cells were clustered based on their transcriptomic profiles (visualized via UMAP ( Figure 3A ) with relative percentages of each cell type depicted in Figure 3B ). The broad cell types were compared via expression of cell-surface marker genes ( Figure 3C ) as well as percentage of spliced mRNA between patient groups ( Figure 3D ). This was followed by careful manual annotation and inspection of canonical marker genes with identification of myeloid cells via differential expression of genes ( Figure 4 ). As explained by Hegde et al. (13), there is substantial plasticity within MDSC subpopulations during sepsis which informs the relationship between MDSCs and terminally differentiated effector cells ( Figure 5A ). Additional marker genes were used to obtain fine-level annotation of myeloid cell types ( Figure 5B ). Importantly, four distinct populations of MDSCs were identified via this approach ( Figure 5C ), three of which were consistent with classically defined E-, PMN-, and M-MDSCs (11).

A novel fourth population was identified in 60% of the late sepsis patients who developed CCI, as well as both of the acutely septic patients who progressed to CCI ( Table 4 , Figure 6A ). This MDSC subpopulation exhibited gene expression patterns that were partially consistent with both M- and PMN-MDSCs. We thus labeled these cells “hybrid” (H)-MDSCs. Although one of the patients with CCI had a much greater number of H-MDSCs than other patients, it should be noted no H-MDSCs were observed in late sepsis patients who had rapidly recovered, or acutely septic patients who progressed to rapid recovery ( Table 4 ). One of our acutely septic patients died within 14 days, giving us an acute sepsis mortality rate of 25%. This is similar to previously reported literature in septic ICU patients (68). The patient that experienced an early death had similar proportions of MDSC subpopulations compared to the other acutely septic patients (M-MDSCs: 61% vs 51 ± 13%, PMN-MDSCs: 0.2% vs 0.7 ± 0.6%, E-MDSCs: 0.08% vs 2 ± 3%, H: 0% vs 0.2 ± 0.5%).

The majority of MDSC-specific genes in septic patients were shared by at least two of the four subpopulations (66%, n=270 genes), with 36% (n=147 genes) significantly expressed by all four ( Figure 6B , Supplementary File 2 ). Although the MDSC subpopulations were fairly similar in terms of overlapping genetic expression, we identified seven genes uniquely expressed by H-MDSCs: RGDP5, TBL1X, MBNL1, SERF2, ATP5F1E, MT-ND1, and MT-ATP6 ( Figure 6C ). Gene expression was downregulated in RANBP2-like and GRIP domain-containing protein 5 (RGDP5), transducin (beta)-like 1X-linked (TBL1X), and muscleblind-like splicing regulator 1 (MBNL1) (69, 70). TBL1X regulates transcriptomic pathways and is upregulated in malignancy (69). MBNL1 regulates alternative splicing and can be up- or downregulated depending on the type of cancer (71). Genes with upregulated expression included small EDRK-rich factor 2 (SERF2), ATP synthase F1 subunit epsilon (ATP5F1E), NADH-ubiquinone oxidoreductase chain 1 (MT-ND1), and mitochondrially encoded ATP synthase membrane subunit 6 (MT-ATP6). The latter three genes encode proteins involved in mitochondrial metabolism and function (72).

We next sought to identify differential genetic expression between our septic cohorts, specifically looking at differences between MDSCs in late sepsis patients who rapidly recovered compared to those who developed CCI. For differential expression across sepsis groups, only M-MDSCs were sufficiently present per group for fitting a linear mixed model with multiple subjects, as MDSCs are a relatively rare population overall. We identified four differentially expressed genes in M-MSDCs using this method: CD163, IER2, CTSZ, and SNX3 ( Figure 6D ). Expression of CD163 (73), a gene responsible for controlling inflammation, was significantly lower in M-MDSCs in late sepsis patients with CCI versus acutely septic patients. SNX3 has been identified as a potential septic biomarker (74), and was significantly upregulated in patients with CCI compared to acutely septic patients. IER2 was significantly higher expressed in late sepsis patients who rapidly recovered compared to acutely septic patients. IER2 is known to be upregulated in response to external stimuli including infection (75, 76). CTSZ expression was significantly higher in patients with CCI compared to patients who rapidly recovered after sepsis, and has been previously identified as a septic marker in mice (77).

The plasticity of the H-MDSC subpopulation is evident in the increased per-cell proportion of unspliced mRNA, indicating more active transcription. Only E-MDSCs had a higher proportion of unspliced mRNA in the myeloid compartment ( Figure 7A ). To examine factors driving cellular activities, we identified genes with a high average proportion of unspliced mRNA within each cell subpopulation and performed enrichment analysis to identify relevant biological processes. Rather than focusing on individual ontologies, we used a network-based approach to cluster similar significantly enriched biological functions for each MDSC subpopulation ( Figures 7B-E ) (52). Not surprisingly, actively transcribed genes in all MDSC subpopulations were enriched for activities pertaining to ‘immune activation.’ While PMN- and M-MDSCs had more biologically distinct functions, H-MDSCs shared enrichment with both cell types, specifically pertaining to pathophysiological septic-related processes including ‘organonitrogen’ and phosphorus-related processes.

Having described the MDSC subpopulations, we next set out to incorporate these findings into differentiation pathways of the myeloid compartment in septic patients. Quantifying transcriptional kinetics via RNA velocity estimation revealed complex, fluid relationships between MDSC phenotypes ( Figure 8A ). As expected, M-MDSCs appeared to serve as the bridge between early immunosuppressive cell types and mature myeloid cells such as monocytes and cDCs ( Figure 8B ). As our analysis was based on PBMCs, it was not possible to compare the transition from MDSCs to mature granulocytes (PMNs). Estimating the graph connectivity between monocyte-lineage cell types allowed us to quantify the strength of each undirected relationship, and showed that MDSC subpopulations are both highly interconnected and much more internally similar to each other than they are to populations of terminally-differentiated myeloid cells ( Figures 8B, C ).

After analyzing velocity-inferred cell state transitions performed with CellRank, all E-, PMN-, and the vast majority of M-MDSCs cell states were classified as progenitor-like or transitioning-like ( Figure 8D ). Only H-MDSCs contained a significant proportion of cells in a plasticity-like state with high probabilities for both initial and terminal cell states (in which cells remain H-MDSCs) ( Figure 8E ) (57). Supporting this, significant variation was observed in the likelihood of an H-MDSC staying an H-MDSC when estimated by absorption probabilities from CellRank, with a mean (SD) probability of 0.33 (0.29) ( Figure 8F ). No other cell types were likely to end up as H-MDSCs. To better characterize the biology underlying commitment to the H-MDSC cell fate, lineage driver genes (genes significantly correlated with the probability of becoming an H-MDSC) were identified by computing Spearman correlations of expression with absorption probabilities. Highly correlated genes were diverse in function and included inflammation-associated genes such as S100A8, -9, and -12, along with immunoregulatory genes ALOX5A, RETN, and IL1R2.

Next, we investigated differences in cell states across sepsis groups for each MDSC subpopulation. As H-MDSCs were not observed in sepsis patients who experienced rapid recovery, they were not included for this analysis. M-MDSCs were highly consistent between septic patients at day 4 and days 14-21 in terms of their cell states and kinetics ( Figure 9A ). Interestingly, PMN-MDSCs displayed the most heterogeneity, specifically in late sepsis patients with CCI compared to both day 4 septic patients and late sepsis patients who rapidly recovered. PMN-MDSCs in late sepsis patients who developed CCI had significantly slower differentiation speed, higher cell state stability, and lower initial state probabilities ( Figure 9B ). This is consistent with PMN-MDSCs persisting in CCI compared to patients who rapidly recover after sepsis. E-MDSCs in late sepsis patients with CCI also showed significantly lower differentiation progression than acutely septic patients or late sepsis patients who rapidly recovered, along with a higher degree of cell commitment along the differentiation trajectory compared to acutely septic patients ( Figure 9C ).

As previously stated, CD66b⁺-isolated PBMCs met the criteria of MDSCs in their ability to suppress either T-lymphocyte cytokine/chemokine production or T-lymphocyte proliferation ex vivo ( Figures 5 , 6 ) (4, 16), although CD66b⁺-isolated PBMCs were not identical in their suppressive activity from acutely septic patients or late time periods after severe infection. Interestingly, whether using cell-surface markers or transcriptomic analysis of the current dataset, differential expression of several key MDSC genes published in the cancer literature did not reach significance and/or were modestly expressed in septic individuals ( Figure 10 ). For example, although there was upregulation of genes in the S100A and MMP superfamilies, differential expression of ARG1, IL-10, NOS2, and TGFB1 did not reach significance (although transcripts from all genes were detected).