Cellular and transcriptional dynamics of human neutrophils at steady state and upon stress

Montaldo, Elisa; Lusito, Eleonora; Bianchessi, Valentina; Caronni, Nicoletta; Scala, Serena; Basso-Ricci, Luca; Cantaffa, Carla; Masserdotti, Alice; Barilaro, Mattia; Barresi, Simona; Genua, Marco; Vittoria, Francesco Maria; Barbiera, Giulia; Lazarevic, Dejan; Messina, Carlo; Xue, Elisabetta; Marktel, Sarah; Tresoldi, Cristina; Milani, Raffaella; Ronchi, Paola; Gattillo, Salvatore; Santoleri, Luca; Di Micco, Raffaella; Ditadi, Andrea; Belfiori, Giulio; Aleotti, Francesca; Naldini, Matteo Maria; Gentner, Bernhard; Gardiman, Elisa; Tamassia, Nicola; Cassatella, Marco Antonio; Hidalgo, Andrés; Kwok, Immanuel; Ng, Lai Guan; Crippa, Stefano; Falconi, Massimo; Pettinella, Francesca; Scapini, Patrizia; Naldini, Luigi; Ciceri, Fabio; Aiuti, Alessandro; Ostuni, Renato

doi:10.1038/s41590-022-01311-1

Download PDF

Resource
Published: 22 September 2022

Cellular and transcriptional dynamics of human neutrophils at steady state and upon stress

Nature Immunology volume 23, pages 1470–1483 (2022)Cite this article

28k Accesses
48 Citations
70 Altmetric
Metrics details

Subjects

Abstract

Traditionally viewed as poorly plastic, neutrophils are now recognized as functionally diverse; however, the extent and determinants of neutrophil heterogeneity in humans remain unclear. We performed a comprehensive immunophenotypic and transcriptome analysis, at a bulk and single-cell level, of neutrophils from healthy donors and patients undergoing stress myelopoiesis upon exposure to growth factors, transplantation of hematopoietic stem cells (HSC-T), development of pancreatic cancer and viral infection. We uncover an extreme diversity of human neutrophils in vivo, reflecting the rates of cell mobilization, differentiation and exposure to environmental signals. Integrated control of developmental and inducible transcriptional programs linked flexible granulopoietic outputs with elicitation of stimulus-specific functional responses. In this context, we detected an acute interferon (IFN) response in the blood of patients receiving HSC-T that was mirrored by marked upregulation of IFN-stimulated genes in neutrophils but not in monocytes. Systematic characterization of human neutrophil plasticity may uncover clinically relevant biomarkers and support the development of diagnostic and therapeutic tools.

Single-cell transcriptome profiling reveals neutrophil heterogeneity in homeostasis and infection

Article 27 July 2020

Xuemei Xie, Qiang Shi, … Hongbo R. Luo

Single-cell profiling reveals the trajectories of natural killer cell differentiation in bone marrow and a stress signature induced by acute myeloid leukemia

Article Open access 25 November 2020

Adeline Crinier, Pierre-Yves Dumas, … Éric Vivier

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

Article 24 April 2023

Andrew J. Kwok, Alice Allcock, … Julian C. Knight

Main

Neutrophils, the most abundant leukocytes in peripheral blood (PB), ensure host immunity by sensing and phagocytosing invading pathogens, as well as by releasing cytotoxic molecules via granule discharge or neutrophil extracellular traps (NETs)¹. Individuals with severe congenital neutropenia or chronic granulomatous disease, diseases characterized by defective neutrophil development or functions, are indeed sensitive to opportunistic infections², and rapid reconstitution of neutrophil counts after hematopoietic stem cell transplantation (HSC-T) is associated with higher survival and hematological recovery in chemotherapy-treated patients³. On the other hand, aberrant neutrophil activation underlies a variety of inflammatory conditions, including autoimmunity, stroke, neurodegeneration and cancer⁴.

The multifaceted activities of neutrophils in health and disease underscore a remarkable functional diversity⁵. In this context, traditional views of neutrophils as short-lived effectors with limited plasticity are challenged by findings that, already at the steady state, these cells persist within organs and acquire tissue-specific genomic programs⁶. Heterogeneity of neutrophils also reflects the output of bone marrow (BM) granulopoiesis, a demand-adapted process sensitive to homeostatic fluctuations⁷, alterations of the hematopoietic niche or changes in the concentration of mediators such as granulocyte colony-stimulating factor (G-CSF)⁸. During stress-induced myelopoiesis, committed precursors and immature neutrophils undergo expansion and premature release in the blood, where they coexist with terminally differentiated subsets^9,10,11. Neutrophil properties are further diversified as cells are exposed to stress-associated stimuli in the circulation or in target organs, leading to the production and release of a spectrum of inflammatory and regulatory products^{12,13,14,15,16}.

A systematic analysis of the phenotypic and transcriptome changes occurring in human neutrophils during inflammation is a prerequisite for the interpretation and rational targeting of these cells’ activities in homeostasis and disease; however, the extent and drivers of neutrophil heterogeneity in humans have remained elusive. We addressed these issues by performing a comprehensive immunophenotype and transcriptome analysis, at the bulk and single-cell level, of human neutrophils and monocytes in healthy controls and in patients undergoing stress-induced myelopoiesis driven by exposure to G-CSF, myeloablative conditioning followed by HSC-T, development of pancreatic ductal adenocarcinoma (PDAC) or infection with severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2).

Results

Dynamics and phenotype of G-CSF-elicited human neutrophils

To characterize stress-elicited neutrophil dynamics in humans, we performed a comprehensive immunophenotype analysis of PB or BM samples from control and age-matched individuals undergoing G-CSF treatment (Supplementary Tables 1–6). Unfractionated blood samples were analyzed using multiparametric flow cytometry with a panel of 16 antibodies able to quantify up to 28 subsets of hematopoietic stem and progenitor cells (HSPCs), committed precursors and differentiated cells¹⁷ (Extended Data Fig. 1a and Supplementary Table 7). To assess neutrophil phenotypes, PB and BM samples were additionally subjected to density gradient separation followed by targeted flow cytometry analyses of CD15⁺CD66b⁺ cells that sedimented in the erythrocyte/granulocyte fraction (referred to as normal-density neutrophils (NDNs)) or in the mononuclear cell layer (low-density neutrophils (LDNs)) (Extended Data Fig. 1b–d and Supplementary Table 7). Neutrophils indeed show variable buoyant densities according to changes in their granule content and nuclear morphology during maturation and/or activation. We observed a robust mobilization of myeloid-biased HSPCs in G-CSF-treated donors (Extended Data Fig. 1e–h), concomitant with a preferential surge of circulating neutrophils and other myeloid cells (Fig. 1a and Extended Data Fig. 1i,j). In line with previous studies^18,19, G-CSF exposure led to increased numbers of neutrophil precursors (SSC^hiCD33/CD66b⁺CD38⁺CD11c⁻CD10⁻) and of immature neutrophils (SSC^hiCD33/CD66b⁺CD38⁻CD11c^−/+CD10⁻) (Fig. 1b,c and Extended Data Fig. 1a,k). G-CSF-treated donors displayed high frequencies of LDNs (Fig. 1d,e) that expressed low levels of the neutrophil differentiation and activation markers CD11b, CD11c, CD62L, CD16 and CD10 (Fig. 1d,f and Extended Data Fig. 2a). G-CSF-elicited LDNs were heterogeneous and included cells corresponding to neutrophil precursors (CD15⁺CD66b⁺CD49d⁺CD16⁻), immature neutrophils (CD15⁺CD66b⁺CD49d⁻CD16^int) and mature neutrophils (CD15⁺CD66b⁺CD49d⁻CD16^hi) (Fig. 1g–i and Extended Data Fig. 2b,c)^9,10,11. Ex vivo assays of 5-ethynyl-2′-deoxyuridine (EdU) incorporation confirmed that LDNs contained neutrophil precursors with proliferative potential^9,10 (Fig. 1j,k and Extended Data Fig. 2d,e). As reported previously¹⁹, NDNs from G-CSF-treated donors were mostly composed of mature neutrophils with an activated phenotype (lower expression of CD10 and CD62L; Fig. 1f and Extended Data Fig. 2a) and increased levels of CD35 (CR1) and CD54 (ICAM-1) (Extended Data Fig. 2f,g) than control NDNs. These data highlight the heterogeneous phenotype of G-CSF-elicited human neutrophils.

**Fig. 1: Dynamics and phenotype of neutrophils elicited by G-CSF.**

Dynamics and phenotype of neutrophils during HSC-T or PDAC

We next profiled neutrophil dynamics in age-matched individuals undergoing emergency myelopoiesis secondary to HSC-T with high-intensity conditioning (Supplementary Tables 1–6). PB and/or BM samples were collected from these patients shortly after treatment (first follow-up, PB collected 16–27 days after HSC-T), at clinical recovery (second follow-up, PB and BM collected 28–40 days after HSC-T) and months after HSC-T (third follow-up, PB and BM collected >180 days after HSC-T) (Extended Data Fig. 2h–j). Flow cytometry analyses of unfractionated PB samples highlighted mobilization of phenotypically defined neutrophil precursors in patients receiving HSC-T at the first or second follow-up (Fig. 2a–c), coinciding with the appearance in the circulation of heterogeneous LDNs (Fig. 2d,e) containing proliferating and non-proliferating precursors as well as immature neutrophils (Fig. 2f,g and Extended Data Fig. 2d,e,k,l). We then analyzed PB samples from treatment-naive patients with locally advanced or metastatic PDAC (n = 19) or with intraductal papillary mucinous neoplasms (IPMN) (n = 15), a type of lesion that often precedes tumor onset²⁰. Total neutrophil counts were largely unaltered in patients with PDAC, whereas those of other hematopoietic cells (HSPCs, monocytes and T lymphocytes) were significantly reduced (Extended Data Fig. 2m–p). Thus, the neutrophil-to-lymphocyte ratio (NLR), a frequently used marker of cancer progression²¹, was higher for patients with PDAC than for healthy controls or patients with IPMN (Extended Data Fig. 2q,r). In line with previous reports²², we observed increased frequencies of circulating LDNs in patients with PDAC (Fig. 2h,i), with this cell population spanning cycling and non-cycling precursors, immature, as well as mature neutrophils (Fig. 2j and Extended Data Fig. 2s–v). We highlight that PDAC patients enrolled in the study are significantly older than healthy donors (Supplementary Tables 1, 3 and 4); however, no increase in LDN frequencies nor NLR values were detected in age-matched patients with IPMN (Fig. 2i and Extended Data Fig. 2r), suggesting that age is not a key determinant of neutrophil dynamics. Our results highlight mobilization of heterogeneous LDNs as a hallmark of stress myelopoiesis induced by G-CSF treatment, HSC-T or PDAC.

Functional properties of G-CSF-elicited neutrophils

We next performed ex vivo experiments to assess reactive oxygen species (ROS) production, neutrophil extracellular trap (NET) release and cytokine synthesis by neutrophil populations isolated from control or G-CSF-treated individuals (Fig. 3a). As compared to those of healthy individuals, NDNs from G-CSF-treated donors displayed a weaker respiratory burst upon treatment with phorbol myristate acetate (PMA) (Fig. 3b,c), possibly reflecting the pre-activated phenotype of these cells in vivo. G-CSF-elicited LDNs also showed lower responses to PMA than controls (Fig. 3b,c), with ROS generation occurring mostly within phenotypically mature cells (Fig. 3d,e). Dose–response experiments revealed that, at a limiting dose of stimulus, LDNs from G-CSF-treated donors were substantially less able to produce ROS than matched NDNs (Fig. 3f). We also found that G-CSF-elicited LDNs underwent PMA-driven NETosis less efficiently than NDNs (Fig. 3g), corroborating the notion that LDNs contain neutrophils that have not fully acquired effector capacities. We next measured the levels of a panel of cytokines, chemokines and growth factors in the culture supernatant of neutrophil populations treated ex vivo with the Toll-like receptor (TLR)-8 agonist resiquimod (R848), a powerful stimulator of cytokine release by human neutrophils²³. A differential biosynthetic capacity emerged when comparing control and G-CSF-elicited NDNs, with the latter cells displaying efficient release of inflammatory molecules upon stimulation (Fig. 3h). Notably, LDNs from G-CSF-treated donors were particularly responsive to TLR ligation and synthesized the highest levels of cytokines such as interleukin (IL)-1β, IL-1RA, G-CSF, CCL2, CCL5 and tumor necrosis factor (TNF)-α (Fig. 3h). These findings suggest that, while lacking at least some effector features of terminally differentiated cells, mobilized immature neutrophils retain immune regulatory capacity via cytokine synthesis and release. Collectively, our data underscore the influence of maturation and exposure to growth factors on the functional effector and the immunoregulatory properties of stress-elicited neutrophil subsets.

**Fig. 3: Functional analysis of G-CSF-elicited neutrophils.**

Bulk transcriptome analysis of monocytes, NDNs and LDNs

To define the gene expression programs of myeloid cells at steady state and after stress, we performed bulk RNA-sequencing (RNA-seq) analyses of NDNs, LDNs and monocytes isolated from the PB of healthy controls (n = 19) and G-CSF-treated donors (n = 17) as well as of patients receiving HSC-T (n = 8) and patients with PDAC (n = 15) or IPMN (n = 14). We also analyzed neutrophil differentiation intermediates from BM samples of controls (n = 3) or patients receiving HSC-T (n = 7), generating a total of 210 RNA-seq samples from 73 individuals (Supplementary Table 8). Cell purity after magnetic bead selection or sorting was consistently higher than 95% (Extended Data Fig. 3a–f). Principal-component analysis (PCA) and unsupervised k-means clustering highlighted clear segregation of monocyte, NDN and LDN transcriptomes (Fig. 4a, Extended Data Fig. 4a–c and Supplementary Table 9). Monocytes were characterized by selective expression of transcripts encoding for known myeloid transcription factors (KLF4, IRF8 and MAFB), scavenger receptors (MARCO and MRC1), components of the antigen presentation machinery (HLA-DMA, HLA-DRA and CD74) and inflammatory cytokines (CCL2 and CXCL10) (Fig. 4a; module 6). Notably, monocytes expressed a gene program (module 5) that was shared with stress-elicited LDNs and that included transcripts encoding for regulators of RNA transcription (POLR1A and POLR2L), translation (EIF2A and EEF2) and ribosome biogenesis (RPL10A, RPS23 and BOP1) (Fig. 4a). LDNs displayed high levels of genes encoding for neutrophil granule proteins (MPO, DEFA4 and ELANE), cell cycle regulators (TOP2A), transcription factors (CEBPE) and surface markers (CEACAM8) (Fig. 4a; modules 2 and 4), and they tended to cluster together with developing neutrophils of the BM from healthy donors (Fig. 4a and Extended Data Fig. 4a–c). LDNs were also characterized by low basal expression of inflammatory response genes (GBP1, OASL, IL1B and TNF) that were instead transcribed in NDNs (Fig. 4a; modules 1 and 3), a finding that was confirmed at the protein level for IL-1β (Extended Data Fig. 5a,b). Collectively, these data indicate that stress-elicited LDNs are characterized by a gene expression program distinct from that of monocytes or NDNs and largely comparable to that of developing neutrophils of the BM.

**Fig. 4: Bulk RNA-seq of neutrophils and monocytes upon G-CSF, HSC-T, IPMN or PDAC.**

Transcriptional responses to stress in NDNs and monocytes

We next set out to define stress-induced transcriptional changes in NDNs and monocytes (due to their heterogeneity, cells corresponding to LDNs were studied at the single-cell level; see below). Analysis of differentially expressed genes (DEGs) (Fig. 4b and Supplementary Tables 10–12) and downstream validation (Extended Data Fig. 5c–g) uncovered a profound transcriptome reprogramming of NDNs from G-CSF-treated donors or patients receiving HSC-T, whereas NDNs from patients with PDAC underwent comparatively small changes. There was a limited overlap between DEGs in NDNs from the various experimental conditions, indicating that transcriptional responses of human neutrophils are largely stress-specific (Fig. 4c and Supplementary Table 13). In line with this notion, exposure to G-CSF led to induction in NDNs of genes belonging to Gene Ontology (GO) categories such as mitochondrial gene translation, oxidative phosphorylation or leukocyte-mediated immunity and to repression of interferon (IFN) responses (Fig. 4d and Supplementary Tables 14–16). On the other hand, NDNs from patients receiving HSC-T showed a clear IFN signature, and they expressed genes of defense response and mitochondrial translation GOs (Fig. 4d and Supplementary Tables 14–16). Increased expression of IFN response genes and of transcripts controlling fatty acid metabolism was also measured in NDNs from patients with PDAC (Fig. 4d and Supplementary Table 14–16). Notably, we found that monocytes from G-CSF-treated donors or patients receiving HSC-T showed limited transcriptional changes as compared to what observed in NDNs from the same individuals (Fig. 4e and Supplementary Tables 17–22). These data highlight a remarkable plasticity of human neutrophils in vivo.

Analysis of plasma factors elicited by G-CSF, HSC-T or PDAC

To identify soluble factors underlying neutrophil dynamics upon stress, we quantified a panel of plasma cytokines, chemokines and growth factors in the PB of subjects enrolled in the study (Supplementary Table 23). G-CSF administration was associated to a drastic upregulation of G-CSF and IL-1RA, as well as to a mild increase of inflammatory cytokines that included IFN-γ, IL-18 and CXCL10 (Fig. 5a,b and Extended Data Fig. 6a). G-CSF treatment also led to lower levels of CXCL12 (Fig. 5a,b), a key BM homing signal for CXCR4⁺ HSPCs and immature myeloid cells. We observed a marked inflammatory skewing of the plasma cytokine profile in patients receiving HSC-T sampled up to 1 month after transplant, with increased levels of factors controlling myeloid cell differentiation (G-CSF, M-CSF and IL-6), recruitment (IL-8, CCL7 and CCL3) and activation (IL-18, IL-12, IL-1α and IL-1β) (Fig. 5a,c). The most upregulated plasma molecules in patients receiving HSC-T were the IFN-stimulated chemokines CXCL9 and CXCL10, in line with a significant elevation of IFN-α2 and IFN-γ shortly after HSC-T (Fig. 5a,c). Patients with PDAC, but not with pre-malignant IPMN, also showed higher levels of proinflammatory cytokines, namely, IL-6, IL-8, CCL3 and M-CSF as well as CXCL9 and CXCL10 (Fig. 5a and Extended Data Fig. 6b). Our data indicate that G-CSF treatment, HSC-T or PDAC development are characterized by a systemic increase in PB concentration of inflammatory molecules known to drive stress myelopoiesis. In this context, plasma levels of G-CSF, IL-6 and IL-8 were positively associated with the frequencies of mobilized neutrophil precursors or LDNs when combining samples from all groups (Fig. 5d,e and Extended Data Fig. 6c–e). A correspondence between plasma cytokine profiles and transcriptional dynamics of neutrophils was evident, as exemplified by the increased levels of IFNs in PB and upregulation of IFN response genes in NDNs from patients receiving HSC-T.

**Fig. 5: Plasma factors underlying stress myelopoiesis upon G-CSF, HSC-T, IPMN or PDAC.**

Transcriptional diversity of human neutrophils upon stress

We next performed single-cell RNA-seq (scRNA-seq) on CD15⁺ cells isolated from PB or BM samples of healthy controls (PB n = 2, BM n = 2), G-CSF-treated donors (PB, n = 4), patients receiving HSC-T (PB n = 3, of which one received G-CSF post-transplant and BM n = 2), or patients with PDAC (PB, n = 5) (Supplementary Table 24). This sorting strategy enabled us to recover the full spectrum of developing neutrophils, from precursors to terminally differentiated cells (Extended Data Fig. 7a). After normalization and filtering, our dataset included transcriptomes from 130,628 cells, of which 1,059 were classified as contaminants (Supplementary Table 25). Graph-based clustering analysis revealed an extensive diversity of human neutrophils, with cells being distributed in the Uniform Manifold Approximation and Projection (UMAP) embedding according to their maturation stage, tissue location, exposure to stress signals and donor/patient identity (Fig. 6a,b, Extended Data Fig. 7b and Supplementary Table 25). We next employed curated gene signatures from developing human neutrophils^24,25 to annotate UMAP clusters (Fig. 6a–c and Extended Data Fig. 7c,d). Cells in cluster 1 expressed the highest levels of a transcriptional module previously associated to neutrophil-committed progenitors²⁵, which include genes encoding for azurophilic granules (MPO and ELANE) and cell cycle proteins (MKI67 and TOP2A); we annotated this cluster as ‘precursors’. We then defined ‘early immature’ (cluster 2 and 3) and ‘immature’ (cluster 4–14) neutrophils based on increasing expression of specific (CAMP, LTF and LCN2) or gelatinase (MMP9 and CTSB) granule genes²⁴, which are progressively transcribed along the transitions from promyelocytes to band cells²⁵. Finally, ‘mature’ neutrophils (cluster 15–24) were defined by expression of SELL, MME and CXCR4 in the BM or NAMPT, CXCR2 and SOD2 in PB (Fig. 6a–c and Extended Data Fig. 7c,d). Mapping of neutrophil gene modules, previously defined by bulk RNA-seq (Fig. 4a), onto single-cell transcriptome data confirmed that LDNs (gene modules 2, 4 and 5) span precursors, early immature, immature and mature neutrophils of the BM, whereas NDNs (module 1) correspond to mature PB neutrophils (Extended Data Fig. 7e). In line with the observed patterns of LDN mobilization, we found that scRNA-seq clusters of precursors, early immature and immature neutrophils were predominantly localized in the BM at steady state but became evident in the PB of G-CSF-treated donors, patients receiving HSC-T or patients with PDAC (Fig. 6a,b and Extended Data Fig. 7f,g). Neutrophils at various developmental stages showed differential patterns of inducible gene expression in response to stress. Mature cells from controls, G-CSF-treated donors, patients receiving HSC-T or patients with PDAC were segregated from each other in scRNA-seq, whereas precursors and immature neutrophils from all experimental conditions tended to cluster together (Fig. 6a,b, Extended Data Fig. 7f,g and Supplementary Table 25). In this context, distinct gene signatures were evident in stress-elicited mature neutrophils. G-CSF exposure was associated with higher expression of transcripts such as SERPINA1, CR1, CX3CR1, CD177, LAIR1 and CD14, whereas neutrophils from a set of patients with PDAC expressed IFN response genes (IRF1 and GBP1) (Fig. 6a–c and Supplementary Table 26). Mature neutrophils from patients receiving HSC-T sampled early after transplant expressed high levels of genes, such as OAS2, CD274, AIM2 and GBP5 (Fig. 6a–c and Supplementary Table 26), which were associated with IFN responses (Fig. 6d). This signature became less evident at later time points (Fig. 6e,f and Supplementary Table 27), consistent with a progressive return to steady state. Collectively, these data indicate that the combined mobilization and exposure to inflammatory factors drive divergent developmental trajectories in stress-elicited neutrophils, resulting in the acquisition of stimulus-specific gene expression programs (for example, IFN signature) in terminally differentiated cells (Extended Data Fig. 7h). More generally, our scRNA-seq analyses uncover a high degree of transcriptional heterogeneity of circulating human neutrophils, dictated by factors such as the differentiation state, their release in circulation and the immunological status of the host.

**Fig. 6: Single-cell RNA-seq analysis of human neutrophils at steady state and upon stress.**

Dynamics of neutrophil differentiation during stress

We next set out to dissect how exposure to stress signals impacted on the continuum of neutrophil differentiation. We first applied cellHarmony²⁶ on single-cell transcriptomes of PB and BM neutrophils from healthy controls to build a reference dataset of steady-state neutropoiesis and to define cell states corresponding to specific developmental intermediates (Supplementary Table 28). Next, we matched scRNA-seq data of neutrophils from PB and, when available, BM samples from G-CSF-treated donors, patients receiving HSC-T and patients with PDAC (termed query) to the previously defined cellHarmony clusters. This approach enables unbiased co-clustering of neutrophils at an equivalent maturation phase and precise quantification of stress-induced transcriptional changes (Fig. 7a,b, Extended Data Fig. 8a–c and Supplementary Tables 29–31). We observed clear differences in both the dynamics and the levels of expression of developmental genes in neutrophils from G-CSF-treated donors and patients receiving HSC-T, with supra-physiological and prolonged expression of marker genes of neutrophil precursors (MPO, DEFA4, ELANE, RNASE2 and PLAC8) (Fig. 7a–c, Extended Data Fig. 8c and Supplementary Tables 29 and 30). Analogously, genes such as TSPO, MMP8, HP, FCN1, FCER1G and S100A6 were hyper-expressed in immature neutrophils from G-CSF-treated donors and patients receiving HSC-T, and they were prematurely and/or persistently transcribed even at earlier or subsequent differentiation stages (Fig. 7a–c, Extended Data Fig. 8c and Supplementary Tables 29 and 30). A similar, although less pronounced, behavior was observed in neutrophils from patients with PDAC (Extended Data Fig. 8a,b and Supplementary Table 31). Our data show that exposure to inflammatory factors leads to substantial changes in the dynamics of expression of neutrophil developmental genes, possibly supporting the enhanced cellular outputs of granulopoiesis during stress.

**Fig. 7: Transcriptome reprogramming of human neutrophils upon stress.**

Transcriptional changes of human neutrophils during stress

To determine how the pre-existing developmental state impacted on stimulus-induced reprogramming of neutrophils, we performed differential gene expression analyses between reference and query datasets in cellHarmony clusters. These studies uncovered sets of genes that were up- or downregulated upon stress in a developmental state-specific manner (Fig. 7d,e, Extended Data Fig. 8d and Supplementary Tables 32–34). In patients receiving HSC-T, transcripts upregulated in differentiated neutrophils (such as ISG15, IFI6 or STAT1) were poorly induced in precursors or in immature neutrophils (Fig. 7d,f and Supplementary Table 33). Conversely, genes induced in precursors were not induced in mature neutrophils (Fig. 7d,f and Supplementary Table 33). The latter behavior was also evident in cells from G-CSF-treated donors or patients with PDAC, with stress-inducible gene expression programs being largely distinct between neutrophils at various developmental states (Fig. 7e,g and Supplementary Table 32). In line with this notion, GOs of cluster-specific genes were distinct. Mature neutrophils from patients receiving HSC-T upregulated genes belonging to IFN response and antiviral defense GOs, whereas precursors and immature cells from the same individuals upregulated genes involved in RNA processing, translation and protein biosynthesis (Fig. 7h and Supplementary Table 35). On the other hand, mature neutrophils from G-CSF-treated donors displayed high expression of genes related to ATP and carbohydrate metabolic process, macrophage activation and cell adhesion (Fig. 7i and Supplementary Table 36). We observed a limited overlap between G-CSF-, HSC-T or PDAC-induced genes in cellHarmony clusters (Extended Data Fig. 8e and Supplementary Table 37), reinforcing the notion that transcriptional responses of neutrophil subsets are stress-specific. We next set out to validate and extend our findings in patients infected with severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), an occurrence associated with emergency granulopoiesis and aberrant neutrophil activation^27,28. Publicly available scRNA-seq datasets of peripheral blood mononuclear cells (PBMCs) from patients with coronavirus disease 19 (COVID-19)²⁷ were integrated with single-cell transcriptomes generated in this study and subjected to cellHarmony and differential gene expression analyses (Extended Data Fig. 9a). In keeping with our previous observations, viral infection altered both the dynamics and the expression levels of developmental genes in neutrophil mobilized upon stress (Extended Data Fig. 9b and Supplementary Table 38). Genes induced in neutrophils from patients with COVID-19 differed for the various developmental intermediates (Extended Data Fig. 9c and Supplementary Table 39), and they were enriched in distinct GO categories (Extended Data Fig. 9d and Supplementary Table 40). Mature neutrophils from SARS-CoV-2-infected individuals displayed an antiviral response signature that was absent or much less evident in precursors or immature neutrophils (Extended Data Fig. 9d). Indeed, differentiated neutrophils from patients with COVID-19 (but not precursors or immature cells) upregulated genes such as IFITM3, LY6E, GBP1, GBP5, IFI6, ISG15 and FECR1G (Extended Data Fig. 9e). Collectively, our data indicate that stress-elicited neutrophils undergo context-dependent transcriptome reprogramming in vivo, in a manner that reflects both the developmental stage and the type of stimuli to which the latter cells are exposed.

Transcriptional responses to IFN by developing neutrophils

Bulk and scRNA-seq analyses highlighted a marked tendency of differentiated human neutrophils to undergo transcriptome reprogramming in response to IFN. These cells showed dynamic expression of IFN-stimulated genes in patients receiving HSC-T, to a degree that even surpassed that of monocytes from the same individuals (Fig. 4d,e and Supplementary Table 18). Furthermore, higher expression of IFN-stimulated genes was observed in mature neutrophils from patients receiving HSC-T (or patients with COVID-19) as compared to less-differentiated cells (Fig. 7b,d,f and Extended Data Fig. 9b–e). To determine how the differentiation stage of neutrophils correlated with their transcriptional responses to IFNs, we performed scRNA-seq on developing neutrophils treated ex vivo with IFN-β or IFN-γ. CD15⁺ cells from cord blood (CB) samples of healthy donors were isolated and pooled to capture the entire spectrum of neutrophil maturation (Extended Data Fig. 10a–e and Supplementary Table 41). We obtained single-cell transcriptomes from 22,440 neutrophils, which clustered in the t-distributed stochastic neighbor embedding (t-SNE) plot according to their developmental stage and the type of treatment (Fig. 8a–d, Extended Data Fig. 10e and Supplementary Table 42). The neutrophil composition of CB largely reflected that of BM and PB, with defined populations of precursors (cluster 1), immature (cluster 2) and mature neutrophils (cluster 3−6) (Fig. 8a,b and Extended Data Fig. 10e). Cluster 6 corresponded to a population of CB neutrophils expressing chemokine genes (CXCL2, CCL3 and CCL4) that was not clearly detectable in BM or PB samples (Extended Data Fig. 10f) and was not investigated further. Treatment of mature neutrophils with IFN-β (cluster 4) or IFN-γ (cluster 5) led to upregulation of antiviral response genes (IFIT1, RSAD2, ISG15 and OASL) or of known IFN-γ target genes (GBP5, CD274, CD69 and SOCS1), respectively (Fig. 8e–g and Supplementary Table 43). On the other hand, precursors (cluster 1) or immature neutrophils (cluster 2) treated with IFNs did not segregate from controls in the scRNA-seq analyses at the clustering resolution used (Fig. 8a–e and Supplementary Table 42). Developing neutrophils showed detectable but weaker induction of IFN-regulated genes than mature cells (Fig. 8f–i), suggesting lower responsiveness to cytokine stimulation. In line with this notion, IFNAR1 and IFNGR1 (transcripts encoding for key IFN-β and IFN-γ receptor subunits, respectively) were highly expressed in mature neutrophils, but not in immature neutrophils or in precursors (Fig. 8j,k). Collectively, our data underscore the extensive transcriptional plasticity of differentiated human neutrophils upon stimulation with environmental agents.

**Fig. 8: Maturation stage-dependent transcriptome dynamics in neutrophils stimulated with IFNs.**

Discussion

In this study, we combined multiparametric immunophenotyping, quantification of plasma factors, bulk and single-cell genomics and computational modeling to dissect cellular dynamics and molecular diversity of human neutrophils at homeostasis and upon stress. Our study extends previous analyses in mice^{14,16,29,30,31}, and it uncovers principles of neutrophil gene expression in relevant conditions of stress-induced myelopoiesis. Treatment with G-CSF of healthy individuals for HSC mobilization, HSC-T in chemotherapy-treated patients and development of pancreatic cancer elicited a common immunological response, namely, release in the blood of developing neutrophils and production of inflammatory cytokines driving neutrophil development and trafficking. The clinical outcomes of the above-described settings are profoundly influenced by the activities of neutrophils, as these cells were shown to control HSC mobilization in response to G-CSF or other agents such as GROβ and AMD3100³², to enable immune protection of the host and vascular repair³³ upon HSC-T and to modulate cancer progression in a context-dependent manner^4,34,35.

The properties of stress-elicited LDNs have been studied in various settings, often by bulk comparison with mature neutrophils³⁶. We report that LDNs are highly heterogeneous and span the entire spectrum of neutrophil differentiation, up to early precursors. In keeping with their immature phenotype, LDNs displayed limited effector properties ex vivo (respiratory burst and NETosis) as compared to terminally differentiated cells. On the other hand, LDNs were particularly efficient at producing cytokines upon TLR ligation, a capacity that was likely supported by the high expression level of transcript and protein biosynthesis genes. The functions of LDNs in vivo remain unclear and include immune modulation or tissue repair³⁷. Our data support the hypothesis that LDNs regulate local and/or systemic inflammation via cytokine production. Furthermore, mobilized LDNs might give rise to mature neutrophils in the periphery and thus support increased cellular demands upon inflammation or damage. In the context of HSC-T, mobilized neutrophil precursors may thus sustain immune reconstitution, and therapeutic approaches that stimulate their production and release could boost recovery from neutropenia following preparative conditioning. Combining lineage tracing with single-cell genomics and functional analyses will elucidate the hierarchy and developmental connections between neutrophil precursors and their progeny, as well as highlight functional implications during homeostasis and disease.

Single-cell transcriptome analysis of mobilized neutrophils exposed to inflammatory stimuli in the blood allowed us to dissect the complex interplay between differentiation and activation in vivo. Stress-elicited neutrophils underwent profound changes in the expression dynamics of developmental genes, and they concomitantly acquired stimulus-specific gene signatures. Notably, both the extent and the type of transcriptional responses to stimulation were different for cells at various maturation stages, with the most evident transcriptome dynamics being observed in differentiated neutrophils. These data support a model whereby combined mobilization and exposure to inflammatory factors elicit divergent neutrophil developmental trajectories that result in the acquisition of context-specific functional programs by terminally differentiated cells. We speculate that an integrated control of developmental and inducible gene expression in neutrophils enables persistent adaptations³⁸, such as those seen in the long-term setting of trained immunity to infection or cancer^39,40.

Transcriptional reprogramming of neutrophils closely mirrored changes of blood cytokine profiles in individuals with stress myelopoiesis. A metabolic and proliferative response underlined neutrophils from G-CSF-treated donors, in line with the known biological actions of the latter molecule. An acute IFN cytokine signature was instead detected in the blood of patients receiving HSC-T early after transplant, possibly reflecting chemotherapy-induced tissue damage, viral reactivation, exposure to pathogens or acute graft versus host disease (GvHD)^41,42; this response was associated with a marked induction of IFN-stimulated genes in circulating neutrophils. Notably, monocytes from the same patients receiving HSC-T underwent minor transcriptional changes, as they upregulated a relatively small set of inflammatory genes. The molecular bases of differential IFN responses by neutrophils and monocytes in vivo remain to be elucidated. A possibility is that lower thresholds of IFN concentrations may be required to drive inducible gene expression in neutrophils; this behavior would be compatible with the existence of differential signal transduction pathways and chromatin dynamics at inflammatory response genes in the two cell types^43,44. We propose that neutrophils act as powerful sensors of environmental stimuli and of IFNs in particular, with the potential to provide accurate transcriptome readouts of signaling networks occurring in the blood. The high responsiveness of neutrophils to IFNs may underlie the relevance of these cells as biomarkers for severe infectious diseases, in line with the reported predictive power of neutrophil gene expression in blood transcriptional signatures of patients with bacterial or viral infection^45,46,47. Future studies will be aimed at determining whether neutrophil transcriptome features can be used as biomarkers of clinical parameters of HSC-T, such as hematopoietic reconstitution, viral reactivation, infections with pathogens or GvHD.

By extending previous efforts to characterize neutrophil properties at the steady state^9,24,25 and in clinically relevant settings, including G-CSF administration¹⁸, lung or heart disease^13,48, viral infection^27,28,49 and cancer^12,15,50, our study represents a step toward a mechanistic understanding of neutrophil diversity in humans. We anticipate that integration of current and future large-scale phenotypic, molecular and functional analyses will enable the development of diagnostic and therapeutic strategies for diseases in which neutrophils are implicated.

Methods

Experimental methods

Study participants and sample collection

Collection of biological samples was compliant with the Declaration of Helsinki and the General Data Protection Regulation, and the study was approved by Ospedale San Raffaele and Azienda Ospedaliera Universitaria Integrata di Verona ethics committee (protocols: TIGET09; MIELO-GEN; NEU-IPMN; and CMRI/55742). A total of 149 participants were enrolled in the study between June 2017 and June 2022. Samples were collected in EDTA-containing sterile vacutainer tubes, stored at 25 °C and processed within 2 hours. Informed consent was obtained by all participants. Participants received no compensation. Age and sex, as well as anonymized clinical information of enrolled participants, are reported in Supplementary Tables 1–6.

Controls and G-CSF-treated donors

Healthy individuals were enrolled at Ospedale San Raffaele and Azienda Ospedaliera Universitaria Integrata di Verona. We collected PB from healthy donors before HSPC mobilization or BM aspiration procedures (n = 55). BM samples were collected from the posterior iliac crests under anesthesia as a standard HSPC donation procedure (n = 14). Mobilized PB was collected from HSPC donors (n = 49) after 5 days of treatment with G-CSF (filgastrim, 10 μg kg⁻¹ per day). CB samples (n = 10) were collected after C-section deliveries at term of gestation of healthy volunteers donating placental tissue.

Patients receiving HSC-T

Patients (n = 16) with hematological malignancies in complete remission were enrolled at Ospedale San Raffaele. They received preparative myeloablative conditioning and underwent a post-transplant pharmacologic prophylaxis regimen to prevent acute and chronic GvHD and infections. Patients underwent allogeneic HSC-T from either haplotype-mismatched related donor (n = 12) or haplotype-matched related donor (n = 4). Fourteen patients received unmanipulated G-CSF-mobilized PB cells and two received unmanipulated BM cells. We collected samples at three time points after HSC-T: first follow-up, early after transplant when white blood cell count reached 500 cells µl⁻¹ for 3 days (PB collected 16–27 days after HSC-T); second follow-up, at clinical recovery (PB and BM collected 28–40 days after HSC-T); and third follow-up, long-term after transplant (PB and BM collected >180 days after HSC-T). Two patients (UPN34 and UPN40) showed delayed or absent engraftment after HSC-T. Among patients receiving post-transplant G-CSF, we only retained UPN47 for scRNA-seq analysis.

Patients with PDAC and IPMN

We collected PB from patients with suspected or proven diagnosis of pre-malignant and malignant lesions of the pancreas at Ospedale San Raffaele. IPMN diagnosis was confirmed by MRI and/or cytological examination on specimens collected via endoscopic ultrasound fine needle aspiration or by histological examination after resection. PDAC diagnosis was confirmed by cytological examination. Samples were retained only for patients with confirmed IPMN (n = 15) or PDAC (n = 19) diagnosis. Exclusion criteria were chemotherapy and/or radiotherapy treatments and occurrence of acute pancreatitis, cholangitis and surgical or invasive endoscopic procedure within 1 month before PB collection.

Cell isolation

Mononuclear cells and granulocytes were separated by density centrifugation over a Lymphoprep (Stemcell Technologies) gradient. PB and CB samples were diluted 1:1 with PBS, and BM and G-CSF-mobilized PB samples were diluted 1:4 with PBS and layered over Lymphoprep. Mononuclear cells were lysed with sterile ACK solution (0.15 M NH₄Cl, 10 mM KHCO₃ and 0.1 mM EDTA) for 5 min at 25 °C to remove residual erythrocytes and counted in the presence of Trypan blue (Sigma) to evaluate cell vitality. Monocytes and LDNs were isolated from the mononuclear cell fraction either by FACS (see below) or by magnetic beads with CD14 Microbeads or CD15 Microbeads (Miltenyi Biotec), respectively. The granulocyte-enriched fraction was further purified over a Hetasep (Stemcell Technologies) gradient followed by erythrocytes lysis and vital count with Trypan blue. NDNs were isolated from total granulocytes by magnetic bead sorting using the Neutrophil Isolation kit (Stemcell Technologies). Alternatively, mononuclear cells and granulocytes were isolated by Ficoll-Paque (GE Healthcare Life Sciences) gradient centrifugation and Dextran (Sigma) gradient, as previously described¹⁹. For Cytochrome C reduction assay and supernatant production, total CD66b⁺ neutrophils were isolated by magnetic bead selection by incubating mononuclear cells or granulocytes with fluorescence-conjugated anti-CD66b monoclonal antibody, followed by incubation with specific anti-fluorochrome microbeads (Miltenyi Biotec) according to the manufacturer’s protocol. Purity of bead-sorted cell subsets was evaluated by flow cytometry analysis. A detailed reagent list is reported in Supplementary Table 44.

Flow cytometry

Whole blood staining

Whole blood flow cytometry analysis was performed as described previously¹⁷. Briefly, 500 µl of PB or 100 µl of BM was incubated with 3 ml or 1 ml, respectively, of ACK solution for 10 min at 25 °C and washed twice with PBS. After a final wash in PBS 1% BSA, cells were resuspended in 100 µl of PBS 1% BSA and incubated with fluorochrome-conjugated antibody mix for 30 min at 25 °C in the dark. Cells were washed, resuspended in 100 µl of PBS 1% BSA and incubated for 15 min in the dark with PI at a final concentration of 0.25 μg ml⁻¹. Samples were acquired at LSR-Fortessa or BD FACSymphony A5 SORP Cytometer (BD Biosciences) using DIVA software v.8.0.2 (BD Biosciences). Data were analyzed using FlowJo software v.10.8.0 (TreeStar). Cell populations were gated as previously described¹⁷ with minor modifications, as reported in Supplementary Table 7 and Extended Data Fig. 1a.

Mononuclear cells and granulocyte staining

Cells were resuspended in PBS containing 1% BSA or 2% FBS and 2 mM EDTA and then incubated with FcR blocking reagent human (Miltenyi Biotec) or with 5% human serum at 25 °C for 5 min. Finally, cells were incubated with fluorochrome-conjugated antibody mix for 20 min at 4 °C in the dark. Cell suspension was washed with PBS 1% BSA and acquired at Navios Flow Cytometer using NAVIOS software v.1.3 (Beckman Coulter), MACSQuant 10 or 16 Analyzers using MACSQuantify software v.2.13 (Miltenyi Biotec). For IL-1β intracellular staining, cells were fixed and permeabilized with IC Fixation Buffer (Thermo Fisher Scientific) and intracellular staining permeabilization buffer (BioLegend) according to manufacturer’s instruction and acquired at FACSCanto II using DIVA software v.8.0.2 (BD Biosciences). Data were analyzed with FlowJo v.10.6.2 (TreeStar).

Fluorescence activated cell sorting

PB monocytes, LDNs and BM developmental intermediates were sorted from the mononuclear cell fraction. Samples were stained as described above and sorted at MoFlo XDP (Beckman Coulter) or FACSAria Fusion (BD Biosciences) cell sorters using Summit software v.5.4 (Beckman Coulter) and DIVA software v.8.0.2 (BD Biosciences), respectively. We sorted monocytes as CD3⁻CD56⁻CD19⁻CD34⁻ (Lin⁻) CD33⁺CD15⁻CD14⁺ cells and LDNs as (Lin⁻) CD33⁺CD14⁻CD15⁺CD193⁻ cells. BM neutrophils were identified as (Lin⁻) CD14⁻CD33⁺CD15⁺CD193⁻ cells and further fractionated into BM1 CD11b⁻CD16⁻ cells, BM2 CD11b⁺CD16⁻ cells, BM3 CD11b⁺CD16^int and BM4 CD11b⁺CD16^hiCD10⁺ cells. For scRNA-seq experiments, neutrophils were isolated from whole blood after lysis with RBC lysis buffer (BioLegend) and sorted as (Lin⁻) CD14⁻CD33⁺CD15⁺ cells. See also Supplementary Table 7. A detailed reagent list is reported in Supplementary Table 44.

EdU incorporation

Mononuclear cells or total granulocytes were plated at 10⁶ cells ml⁻¹ with RPMI + 10% FBS + 1% Gln + 1% pen/strep in the absence or in the presence of 10 µM EdU. After 18 hours of culture, cells were collected, washed with PBS + 1% BSA, incubated with Fc blocking reagents (Miltenyi Biotec) and stained. Cells were fixed, permeabilized and incubated with reaction cocktail according to Click-iT Plus EdU Flow Cytometry Assay kit (Thermo Fisher Scientific). Samples were acquired at Navios Flow Cytometer using NAVIOS software v.1.3 (Beckman Coulter) and analyzed with FlowJo v.10.6.2 (TreeStar). A detailed reagent list is reported in Supplementary Table 44.

ROS production

Cytochrome C reduction assays or neutrophil/monocyte respiratory burst assay kits (Cayman Chemical) were used. Freshly isolated CD66b⁺ LDNs and/or NDNs were washed and resuspended at 2 × 10⁶ cells ml⁻¹ in Hank’s Balanced Salt Solution (HBSS), pH 7.4, supplemented with 10% FBS, 0.5 mM CaCl₂ and 1 mg ml⁻¹ glucose. O₂⁻ production in response to 20 ng ml⁻¹ PMA (Sigma) was assessed by the Cytochrome C reduction assay (Cayman), as previously described⁵¹. For flow cytometry analysis of ROS, 1 × 10⁵ mononuclear cells or granulocytes were incubated with dihydrorhodamine-123 (Cayman Chemical) and left untreated or stimulated with PMA 20 ng ml⁻¹ for 15 or 30 min. Cells were stained with fluorochrome-conjugated antibodies as described above, acquired at FACSCanto II using DIVA software v.8.0.2 (BD Biosciences) and analyzed with FlowJo v.10.6.2 (TreeStar). LDNs and NDNs were identified after gating on Lin⁻CD15⁺ cells in the PBMC and granulocyte fraction, respectively. A detailed reagent list is reported in Supplementary Table 44.

NETosis

The NETosis assay kit (Cayman) was used. Bead-sorted NDNs and LDNs were resuspended at 1 × 10⁶ cells ml⁻¹ and left untreated or stimulated with PMA 20 nM and incubated at 37 °C for 2 hours. Culture supernatants were removed, and wells were washed to remove soluble elastase. After treatment with S7 nuclease to induce the release of NET-associated elastase, supernatants were collected and elastase activity was evaluated according to manufacturer’s instructions. A detailed reagent list is reported in Supplementary Table 44.

Ex vivo stimulation of NDNs and LDNs

Purified LDNs and NDNs were plated at 5 × 10⁶ ml⁻¹ in the presence of RPMI 1640 medium supplemented with 10% FBS and treated or not with 5 μM R848 (InvivoGen). After 20 hours of culture at 37 °C, neutrophils were collected and spun at 300 g for 5 min. Cell-free supernatants were immediately frozen and stored at −80 °C until use. A detailed reagent list is reported in Supplementary Table 44.

Plasma collection

An aliquot of 300 µl of blood collected into EDTA tubes was centrifuged 5 min at 10,000 g. Plasma was transferred into a clean tube and re-centrifuged 5 min at 10,000 g. Plasma was frozen and stored at −20 °C until use.

ELISA

Cytokine and chemokine concentration in culture supernatants or plasma was measured using Bio-Plex Pro Human Cytokine Screening Panel, 48-Plex (Bio-Rad) according to the manufacturer’s instructions. Acquisition was performed using Luminex instruments and analyzed with Bio-plex manager (Bio-Rad) software. A detailed reagent list is reported in Supplementary Table 44.

Cytospin and May-Grünwald Giemsa staining

We resuspended 100,000 cells in 200 ml of PBS + 2% FBS and deposited on a slide with a Cytospin 4 centrifuge (Thermo Fisher Scientific). Slides were dried for 30 min at 25 °C and stained with May-Grünwald solution (Carlo Erba) for 5 min. After washing with water, slides were stained with Giemsa (Merck) working solution (Giemsa solution diluted 1:10) for 15 min and washed with water. Slides were dried in upright position at 25 °C. Images were acquired in bright field using an Eclipse (Nikon) microscope and NIS-Elements 4.0 software. A detailed reagent list is reported in Supplementary Table 44.

Real-time quantitative PCR

RNA was extracted using the ReliaPrep RNA Cell Miniprep System (Promega) and measured with Qubit RNA HS Assay kit using a Qubit 3.0. Then, 0.5 ng of RNA were retrotranscribed with SuperScript II and cDNA was PCR-amplified with KAPA HiFi HotStart. Target genes amplification was performed with Fast SYBR Green Master Mix on a ViiA 7 Real-Time PCR System. A detailed reagent list is reported in Supplementary Table 44.

Ex vivo stimulation of CB neutrophils

We isolated mononuclear cells and granulocytes, as reported above, from three different CB samples. From the mononuclear cell fraction, we isolated LDNs by performing a double round of magnetic bead sorting using a Neutrophil Isolation kit (Stemcell Technologies). From the granulocyte cell fraction, we isolated NDNs by performing a single round of magnetic bead sorting using Neutrophil Isolation kit (Stemcell Technologies). To ensure a sufficient representation of neutrophil precursors (less abundant cell population) and of NDNs (less efficiently detected by droplet-based scRNA-seq due to low RNA content), LDNs and NDNs from each CB sample were mixed in a ratio of 1:3. We plated a LDN–NDN mix at 10⁶ cells ml⁻¹ in RPMI 1640 + 10% FBS + 1% Gln + 1% pen/strep alone or with G-CSF, IFN-β or IFN-γ all used at 10 ng ml⁻¹. After 4 hours, cells were collected, washed and counted, and for each condition, we mixed cells from different CB in a ratio 1:1:1. The pooled samples were processed for scRNA-seq as described below. A detailed reagent list is reported in Supplementary Table 44.

Bulk RNA sequencing

We extracted total RNA using the ReliaPrep RNA Cell Miniprep System (Promega). RNA concentration was measured with Qubit RNA HS Assay kit using Qubit 3.0, and RNA integrity was evaluated with Agilent RNA 6000 Pico kit using Bionalyzer (Agilent). RNA-seq libraries were generated using the Smart-seq2 method⁵² starting from 0.5 ng of RNA. Retro-transcription was performed using SuperScript II Reverse Transcriptase, and complementary DNA was PCR-amplified (18 cycles) with KAPA HiFi HotStart and purified with AMPure XP beads. After purification, we determined cDNA concentration using Qubit dsDNA HS Assay kit at Qubit 3.0, and we assessed size distribution at Agilent 4200 TapeStation system. We performed the tagmentation reaction starting from 0.5 ng of cDNA for 30 min at 55 °C, and we performed enrichment PCR using 12 cycles. Libraries were purified with AMPure XP beads, quantified using Qubit 3.0 and assessed for fragment size distribution on an Agilent 4200 TapeStation system. Libraries were sequenced on an Illumina NextSeq500 or NovaSeq6000 (single-end, 75-bp read length) according to the manufacturer’s instructions. A detailed reagent list is reported in Supplementary Table 44.

Single-cell RNA sequencing

We isolated total CD15⁺ cells and LDNs (from one G-CSF-stimulated donor) by cell sorting. We generated scRNA-seq libraries using the microfluidics-based approach of Chromium Single-Cell Controller (10x Genomics) using the Chromium Single Cell 3′ Reagent kit v.3.0 according to the manufacturers’ instructions. In each experiment, we loaded sample to obtain a target cell recovery of 10,000 cells. cDNA amplification was performed using 13 PCR cycles. The concentration of the scRNA-seq libraries was determined using Qubit dsDNA HS Assay kit at Qubit 3.0, and size distribution was assessed using an Agilent 4200 TapeStation system. Libraries were sequenced on an Illumina NextSeq500 or NovaSeq6000 instruments (paired-end, 150-bp read length) according to the manufacturer’s instruction. A detailed reagent list is reported in Supplementary Table 44.

Computational methods

Bulk RNA-seq analyses on NDNs, LDNs and monocytes

Data processing

Single-end reads (75 bp) were mapped to the GRCh38 reference genome using STAR aligner (v.2.6.0a)⁵³. The FeatureCounts function from Rsubread package (v.3.7)⁵⁴ was then used to summarize the aligned reads to NCBI Homo sapiens RefSeq genes (hg38) while setting the minMQS option to 3. Downstream analyses on the count matrix of expressed genes (25,064 genes and 210 samples) were performed in R environment (v.4.0.1) with edgeR R package (v.3.20.7)⁵⁵. First, genes with more than one count-per-million (CPM) in at least 15% of the total set of samples (NDNs, LDNs and monocytes) were retained for a total of 8,419 genes and 210 samples. Read counts of expressed genes were then normalized with the trimmed mean of m-values method⁵⁶ using calcNormFactors function. The weighted likelihood empirical Bayes method^57,58 was used to calculate the posterior dispersion estimates through the estimateDisp function. The ComBat_seq function⁵⁹ from the sva package (v.3.38.0)⁶⁰ was used to model and correct the batch effects between the sequencing runs. The PCA of the samples was performed based on the batch-corrected reads per count million.

Heat map of variable genes

Log₂ (CPM + 1) were calculated from the batch-corrected counts and used to compute the gene-wise variance across all samples. The values above the 80th percentile of the resulting variance distribution were selected and the corresponding genes used to perform the unsupervised k-means cluster analysis on the standardized expression values with k equal to six. Hierarchical cluster analysis was then performed on the gene modules and samples using the Pearson correlation as distance method and the ward.D2 agglomerative algorithm as hierarchical clustering method.

Differential gene expression analysis

LDNs (82 samples) were removed from the raw count matrix generated with the FeatureCounts R function. The resulting matrix was composed by 25,064 genes and 128 samples: 70 NDNs and 58 monocytes. Genes with more than one CPM in at least 15% of the NDNs or monocytes were selected. The resulting matrix was composed by 8,362 genes and 128 samples. Read counts were normalized and corrected for batch effects as above. Differential gene expression analysis of myeloid cells after stress with respect to the steady state was performed with edgeR for NDNs and monocytes independently starting from the adjusted count matrix containing both cell types. NDNs and monocytes were selected, and the two datasets were further divided into three stress-related/steady-state datasets, each composed by samples from one of the stress condition (G-CSF, HSC-T and PDAC) and samples at steady state. Only genes with CPM > 1 in at least 30% of the samples composing each sub-dataset were retained. The differential gene expression analysis for each stress and cell type was performed by fitting a negative binomial generalized linear model with robust hyperparameter estimation^57,61 using the glmQLFit function and after computing the dispersion with estimateDisp function. A quasi-likelihood F-test^62,63 was then performed using the glmQLFTest function. The sequencing run ID was included in the design matrix of each comparison as a covariate. Genes with abs(log₂FC) ≥ 1.5 and FDR < 0.05 were considered to be differentially expressed.

Principal-component analysis

PCA of NDNs, LDNs and monocytes was performed on expressed genes with CPM > 1 in at least 30% of the total samples of each cell type and with a variance greater than the 95th percentile of the distribution of gene-wise computed variances.

Gene set enrichment analyses

For each stress condition and cell type, DEGs were ranked by decreasing order of log₂FC in stress versus steady state. Gene set enrichment analyses (GSEA) (v.4.0.3)⁶⁴ was performed on ranked gene lists using GO Biological Process Ontology (c5.go.bp.v.7.4) as gene sets, with number of permutations equal to 1,000.

Single-cell RNA-seq analyses of PB and BM neutrophils

Data processing

Fastq files were generated from raw Illumina BCL files using CellRanger v.6.0.2 (10x Genomics) with CellRanger mkfastq and default parameters. CellRanger count was then used to align sequencing reads to the reference transcriptome GRCh38, to perform unique molecular identifier (UMI) filtering and barcode and UMI counting. Only confidently mapped reads with valid barcodes, UMIs and non-PCR duplicated were retained by the tool. The overall sequencing quality was evaluated by looking at the summary metrics of the web_summary.html file generated by the CellRanger pipeline for each sample. The Seurat v4.0.5R package (https://satijalab.org/seurat/) was then used to perform all downstream analysis. First, we removed cells expressing fewer than 300 unique genes and genes expressed in fewer than three cells from the non-normalized UMI count matrix of each sample. Raw count matrices of all samples were then combined in a single Seurat object (17,625 genes and 143,485 cells) with the use of the merge function. A cell/gene quality control was then performed. We jointly examined the distribution of the count depth (number of counts per barcode) of the number of genes per barcode and of the fraction of counts from mitochondrial genes per barcode. Outlier peaks were then filtered out by thresholding. Cells with a total number of detected molecules <500, indicating low-quality cells or empty droplets, were discarded. We also removed cells with a percentage of reads that map to the mitochondrial genes greater than 10% and cells with a number of detected genes >4,000. The two filters were respectively used to remove low-quality/dying cells and cell doublets or multiplets with an aberrantly high gene count⁶⁵. We also applied a gene-wise filter on the average counts to remove low-abundance genes⁶². The filter threshold was established looking at the distribution of the average counts. Genes with a value less than the 15th percentile of the distribution were removed. The final raw count matrix was composed of 15,020 genes and 130,628 cells. We then applied the sctransform normalization⁶⁶ (SCTransform function) while adjusting for the following confounding sources of variation: the mitochondrial mapping percentage and the cell cycle scores computed with the CellCycleScoring function. Data were then scaled with ScaleData, and the top 1,000 variable features were selected with the ‘vst’ method of the FindVariableFeatures function. A shared nearest neighbor graph was constructed using the FindNeighbors function taking as input the first 50 principal components, computed with the RunPCA function. Cell clusters were defined using a resolution of 1.5, were calculated with the FindCluster function and were visualized in two dimensions using UMAP⁶⁷. Cluster-specific marker genes were identified using the MAST method⁶⁸ through the FindMarkers function with option only.pos = TRUE, min.pct = 0.1 and setting a cutoff of FDR < 0.05.

The scRNA-seq data of patients with COVID-19 generated by Schulte-Schrepping²⁷ were downloaded from FASTGenomics database at https://beta.fastgenomics.org/datasets/detail-dataset-7656cfe94fb14a01b787f4774e555036. The dataset used in our analysis was PBMC 10x from cohort2 (Bonn cohort) composed of 46,611 genes and 3,154 cells relative to 22 patients with COVID-19. From the pre-analyzed seurat_COVID19_Neutrophils_cohort1_10x_jonas_FG_2020-08-19.rds file, we extracted the raw counts and re-analyzed the data by applying the quality control criteria used for our datasets to ensure methodological consistency however conditioned to the distribution and features of the data. We first removed the cells expressing fewer than 300 unique genes and genes expressed in fewer than three cells from the non-normalized UMI count matrix, resulting in 13,957 genes across 3,138 cells. Based on the visual inspection of the distribution of the detected molecules across the retained cells, we removed cells with fewer than 500 detected transcripts indicating low-quality cells or empty droplets. We also removed cells with more than 10% mitochondrial reads and with >2,000 detected genes, indicating putative doublets or multiplets. Genes with few counts (fewer than the 15th percentile based on the distribution of the average gene-wise counts across all cells) were considered uninformative and removed. According to the applied criteria for the quality control of cells and genes, the dataset was finally composed of 12,113 genes and 2,990 cells. On these data, we performed the normalization, the identification of the highly variable features, the scaling, the linear dimensionality reduction and the clustering as described above. A batch-effect correction on the normalized expression matrix was performed to run cellHarmony, using ComBat from the sva package to adjust for potential batch effects between donors.

Gene set enrichment analysis

The top 50 marker genes were ranked by decreasing order of log₂FC > 0. GSEA (v.4.0.3)⁶⁴ was performed on ranked gene lists using GO Biological Process Ontology (c5.go.bp.v.7.4) as gene sets, with number of permutations equal to 1,000.

cellHarmony analyses

scRNA-seq raw count matrices of G-CSF-treated donors, patients receiving HSC-T, patients with PDAC and PB or BM healthy donors (HDs) were merged for each condition and preprocessed and normalized with Seurat v.4.0.5 using the same criteria and methods as described above with the following exceptions: cells with a percent of mitochondrial genes greater than 25%, 10%, 15% and 10% relative to G-CSF, HSC-T, PDAC and HDs, respectively, were removed. The threshold for putative doublets and multiplets was also changed and established to be 3,500 for G-CSF and 4,500 for PDAC after the joint visualization of the number of genes and counts. It remained unchanged for the HSC-T and HD datasets. A batch-effect correction was additionally applied to the normalized count matrix of each dataset using ComBat^69,70 from the sva package (v.3.38.0) to adjust for potential batch effects between donors of the same condition. cellHarmony²⁶ was then applied to match cells at the same differentiation stage between the healthy condition (the reference) and the stress (the query). First, the reference dataset (15,851 HD cells: 10,173 BM cells and 5,678 PB cells) was subjected to an unsupervised analysis with ICGS v.2 (AltAnalyze v.2.1.2) that identified eight distinct clusters corresponding to discrete differentiation stages of BM and blood neutrophils. Two of them were considered contaminants and removed. Options were accepted by default except for the number of ICGS cluster (k) that was set to 15 and the column clustering method that was ‘hopach’. Cells from each stress condition (G-CSF, 30,787 cells; HSC-T, 39,479 cells; PDAC, 21,153 cells; and COVID-19, 2,990 cells) were then matched to the reference with cellHarmony to identify analogous differentiation stages. Pairwise differential gene expression analysis between the query cells and the reference cells was performed for each cluster and for each stress independently with FindMarkers function of Seurat v.4.0.5 R package using the MAST method on jointly preprocessed and SCT-normalized expression matrices (steady state + G-CSF; steady state + HSC-T; steady state + PDAC; and steady state + COVID-19). The minimum detection rate (min.pct) was set to 20%. Genes with log₂FC ≥ 1 and FDR < 0.05 were further considered to be differentially expressed.

Gene set enrichment analysis

Due to the small gene set size of the gene lists generated by applying the log₂FC ≥ 1 threshold, the full-length gene lists previously identified with FindMarkers by applying only the detection rate cutoff of 20% were used to run the GSEA. Genes were ranked by decreasing order of log₂FC in stress versus healthy for each cluster of differentiation. GSEA was performed on ranked gene lists using GO Biological Process Ontology (c5.go.bp.v.7.4) as gene sets, with number of permutations equal to 1,000.

Single-cell RNA-seq analyses of CB neutrophils

Chromium scRNA-seq raw data were preprocessed with CellRanger v.6.0.2 (10x Genomics) as described above. Filtered UMI count matrices of CB neutrophils unstimulated (control), stimulated with IFN-β, IFN-γ and G-CSF were analyzed with Seurat v.4.0.5 R package. Data were first subjected to quality control and cells, and genes were selected/removed based on the same criteria described above (min.cells = 3; min.features = 300; percent.MT < 10; nFeature_RNA < 4,000; and nCount_RNA > 500). The 20th percentile of the overall distribution of gene expression levels was used as threshold to remove poorly expressed genes. Data (13,813 genes and 22,440 cells) were then SCT-normalized and scaled while adjusting for cell-cycle effects and the mitochondrial percentage. The top 1,000 variable features were selected with the ‘vst’ method and used as input for PCA. A shared nearest neighbor graph was constructed using the FindNeighbors function taking as input the first 50 principal components, computed with RunPCA function. Cell clusters were defined using a resolution of 0.3, calculated with the FindCluster function and were visualized in two dimensions using t-SNE. Cluster-specific marker genes were identified using the MAST method through the FindMarkers function. Only genes expressed in at least 10% of either of the two groups were tested.

Statistics and reproducibility

No statistical method was used to predetermine sample size. No data were excluded from the analysis. Datasets used for the specific analyses are reported in Methods. Statistical assumptions, including data distribution, independence of observations and homogeneity of variance, were considered for each dataset, and statistical tests were performed accordingly. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.

Reporting summary

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

Data availability

Bulk and scRNA-seq data generated in this study have been deposited in ArrayExpress under accession nos. E-MTAB-11190 and E-MTAB-11188. scRNA-seq data of patients with COVID-19 ²⁷ were downloaded from FASTGenomics database at https://beta.fastgenomics.org/datasets/detail-dataset-7656cfe94fb14a01b787f4774e555036. Source data are provided with this paper.

References

Ley, K. et al. Neutrophils: new insights and open questions. Sci. Immunol. https://doi.org/10.1126/sciimmunol.aat4579 (2018).
Skokowa, J., Dale, D. C., Touw, I. P., Zeidler, C. & Welte, K. Severe congenital neutropenias. Nat. Rev. Dis. Prim. 3, 17032 (2017).
Article PubMed Google Scholar
Copelan, E. A. Hematopoietic stem-cell transplantation. N. Engl. J. Med. 354, 1813–1826 (2006).
Article CAS PubMed Google Scholar
Jaillon, S. et al. Neutrophil diversity and plasticity in tumour progression and therapy. Nat. Rev. Cancer 20, 485–503 (2020).
Article CAS PubMed Google Scholar
Ng, L. G., Ostuni, R. & Hidalgo, A. Heterogeneity of neutrophils. Nat. Rev. Immunol. 19, 255–265 (2019).
Article CAS PubMed Google Scholar
Ballesteros, I. et al. Co-option of neutrophil fates by tissue environments. Cell 183, 1282–1297 (2020).
Article CAS PubMed Google Scholar
Casanova-Acebes, M. et al. Rhythmic modulation of the hematopoietic niche through neutrophil clearance. Cell 153, 1025–1035 (2013).
Article CAS PubMed PubMed Central Google Scholar
Manz, M. G. & Boettcher, S. Emergency granulopoiesis. Nat. Rev. Immunol. 14, 302–314 (2014).
Article CAS PubMed Google Scholar
Kwok, I. et al. Combinatorial single-cell analyses of granulocyte-monocyte progenitor heterogeneity reveals an early uni-potent neutrophil progenitor. Immunity 53, 303–318 (2020).
Article CAS PubMed Google Scholar
Evrard, M. et al. Developmental analysis of bone marrow neutrophils reveals populations specialized in expansion, trafficking, and effector functions. Immunity 48, 364–379 (2018).
Article CAS PubMed Google Scholar
Dinh, H. Q. et al. Coexpression of CD71 and CD117 Identifies an early unipotent neutrophil progenitor population in human bone marrow. Immunity 53, 319–334 (2020).
Article CAS PubMed PubMed Central Google Scholar
Zilionis, R. et al. Single-cell transcriptomics of human and mouse lung cancers reveals conserved myeloid populations across individuals and species. Immunity 50, 1317–1334 (2019).
Article CAS PubMed PubMed Central Google Scholar
Calcagno, D. M. et al. The myeloid type I interferon response to myocardial infarction begins in bone marrow and is regulated by Nrf2-activated macrophages. Sci. Immunol. https://doi.org/10.1126/sciimmunol.aaz1974 (2020).
Xie, X. et al. Single-cell transcriptome profiling reveals neutrophil heterogeneity in homeostasis and infection. Nat. Immunol. 21, 1119–1133 (2020).
Article CAS PubMed PubMed Central Google Scholar
Zhu, Y. P. et al. CyTOF mass cytometry reveals phenotypically distinct human blood neutrophil populations differentially correlated with melanoma stage. J. Immunother. Cancer https://doi.org/10.1136/jitc-2019-000473 (2020).
Khoyratty, T. E. et al. Distinct transcription factor networks control neutrophil-driven inflammation. Nat. Immunol. 22, 1093–1106 (2021).
Article CAS PubMed PubMed Central Google Scholar
Basso-Ricci, L. et al. Multiparametric whole blood dissection: a one-shot comprehensive picture of the human hematopoietic system. Cytom. A 91, 952–965 (2017).
Article CAS Google Scholar
Pedersen, C. C. et al. Changes in gene expression during G-CSF-induced emergency granulopoiesis in humans. J. Immunol. 197, 1989–1999 (2016).
Article CAS PubMed Google Scholar
Marini, O. et al. Mature CD10(+) and immature CD10(−) neutrophils present in G-CSF-treated donors display opposite effects on T cells. Blood 129, 1343–1356 (2017).
Article CAS PubMed Google Scholar
Crippa, S. et al. Low progression of intraductal papillary mucinous neoplasms with worrisome features and high-risk stigmata undergoing non-operative management: a mid-term follow-up analysis. Gut 66, 495–506 (2017).
Article PubMed Google Scholar
Howard, R., Kanetsky, P. A. & Egan, K. M. Exploring the prognostic value of the neutrophil-to-lymphocyte ratio in cancer. Sci. Rep. 9, 19673 (2019).
Article CAS PubMed PubMed Central Google Scholar
Zhang, J. et al. CD13(hi) neutrophil-like myeloid-derived suppressor cells exert immune suppression through arginase 1 expression in pancreatic ductal adenocarcinoma. Oncoimmunology 6, e1258504 (2017).
Article PubMed PubMed Central Google Scholar
Tamassia, N. et al. Induction of OCT2 contributes to regulate the gene expression program in human neutrophils activated via TLR8. Cell Rep. 35, 109143 (2021).
Article CAS PubMed Google Scholar
Grassi, L. et al. Dynamics of transcription regulation in human bone marrow myeloid differentiation to mature blood neutrophils. Cell Rep. 24, 2784–2794 (2018).
Article CAS PubMed PubMed Central Google Scholar
Calzetti, F. et al. CD66b(−)CD64(dim)CD115(−) cells in the human bone marrow represent neutrophil-committed progenitors. Nat. Immunol. 23, 679–691 (2022).
Article CAS PubMed Google Scholar
DePasquale, E. A. K. et al. cellHarmony: cell-level matching and holistic comparison of single-cell transcriptomes. Nucleic Acids Res. 47, e138 (2019).
Article CAS PubMed PubMed Central Google Scholar
Schulte-Schrepping, J. et al. Severe COVID-19 is marked by a dysregulated myeloid cell compartment. Cell 182, 1419–1440 (2020).
Article CAS PubMed PubMed Central Google Scholar
Silvin, A. et al. Elevated calprotectin and abnormal myeloid cell subsets discriminate severe from mild COVID-19. Cell 182, 1401–1418 (2020).
Article CAS PubMed PubMed Central Google Scholar
Grieshaber-Bouyer, R. et al. The neutrotime transcriptional signature defines a single continuum of neutrophils across biological compartments. Nat. Commun. 12, 2856 (2021).
Article CAS PubMed PubMed Central Google Scholar
Muench, D. E. et al. Mouse models of neutropenia reveal progenitor-stage-specific defects. Nature 582, 109–114 (2020).
Article CAS PubMed PubMed Central Google Scholar
Alshetaiwi, H. et al. Defining the emergence of myeloid-derived suppressor cells in breast cancer using single-cell transcriptomics. Sci. Immunol. https://doi.org/10.1126/sciimmunol.aay6017 (2020).
Hoggatt, J. et al. Rapid mobilization reveals a highly engraftable hematopoietic stem cell. Cell 172, 191–204 (2018).
Article CAS PubMed Google Scholar
Bowers, E. et al. Granulocyte-derived TNF-α promotes vascular and hematopoietic regeneration in the bone marrow. Nat. Med. 24, 95–102 (2018).
Article CAS PubMed Google Scholar
Shaul, M. E. & Fridlender, Z. G. Tumour-associated neutrophils in patients with cancer. Nat. Rev. Clin. Oncol. 16, 601–620 (2019).
Article PubMed Google Scholar
Quail, D. F. et al. Neutrophil phenotypes and functions in cancer: a consensus statement. J. Exp. Med. https://doi.org/10.1084/jem.20220011 (2022).
Cassatella, M. A. & Scapini, P. On the improper use of the term high-density neutrophils. Trends Immunol. 41, 1059–1061 (2020).
Article CAS PubMed Google Scholar
Veglia, F., Sanseviero, E. & Gabrilovich, D. I. Myeloid-derived suppressor cells in the era of increasing myeloid cell diversity. Nat. Rev. Immunol. 21, 485–498 (2021).
Article CAS PubMed PubMed Central Google Scholar
Natoli, G. & Ostuni, R. Adaptation and memory in immune responses. Nat. Immunol. 20, 783–792 (2019).
Article CAS PubMed Google Scholar
Kalafati, L. et al. Innate immune training of granulopoiesis promotes anti-tumor activity. Cell 183, 771–785 (2020).
Article CAS PubMed PubMed Central Google Scholar
Moorlag, S. et al. BCG vaccination induces long-term functional reprogramming of human neutrophils. Cell Rep. 33, 108387 (2020).
Article CAS PubMed PubMed Central Google Scholar
Hill, G. R. & Koyama, M. Cytokines and costimulation in acute graft-versus-host disease. Blood 136, 418–428 (2020).
Article PubMed PubMed Central Google Scholar
Tomblyn, M. et al. Guidelines for preventing infectious complications among hematopoietic cell transplantation recipients: a global perspective. Biol. Blood Marrow Transpl. 15, 1143–1238 (2009).
Article CAS Google Scholar
Tamassia, N. et al. Cutting edge: an inactive chromatin configuration at the IL-10 locus in human neutrophils. J. Immunol. 190, 1921–1925 (2013).
Article CAS PubMed Google Scholar
Tamassia, N. et al. The MyD88-independent pathway is not mobilized in human neutrophils stimulated via TLR4. J. Immunol. 178, 7344–7356 (2007).
Article CAS PubMed Google Scholar
Berry, M. P. et al. An interferon-inducible neutrophil-driven blood transcriptional signature in human tuberculosis. Nature 466, 973–977 (2010).
Article CAS PubMed PubMed Central Google Scholar
Dunning, J. et al. Progression of whole-blood transcriptional signatures from interferon-induced to neutrophil-associated patterns in severe influenza. Nat. Immunol. 19, 625–635 (2018).
Article CAS PubMed PubMed Central Google Scholar
Aschenbrenner, A. C. et al. Disease severity-specific neutrophil signatures in blood transcriptomes stratify COVID-19 patients. Genome Med. 13, 7 (2021).
Article CAS PubMed PubMed Central Google Scholar
Song, D. et al. A cellular census of human peripheral immune cells identifies novel cell states in lung diseases. Clin. Transl. Med. 11, e579 (2021).
Article CAS PubMed PubMed Central Google Scholar
Bost, P. et al. Deciphering the state of immune silence in fatal COVID-19 patients. Nat. Commun. 12, 1428 (2021).
Article CAS PubMed PubMed Central Google Scholar
Hsu, B. E. et al. Immature low-density neutrophils exhibit metabolic flexibility that facilitates breast cancer liver metastasis. Cell Rep. 27, 3902–3915 (2019).
Article CAS PubMed Google Scholar
Serra, M. C., Calzetti, F., Ceska, M. & Cassatella, M. A. Effect of substance P on superoxide anion and IL-8 production by human PMNL. Immunology 82, 63–69 (1994).
CAS PubMed PubMed Central Google Scholar
Picelli, S. et al. Full-length RNA-seq from single cells using Smart-seq2. Nat. Protoc. 9, 171–181 (2014).
Article CAS PubMed Google Scholar
Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
Article CAS PubMed Google Scholar
Liao, Y., Smyth, G. K. & Shi, W. The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res. 47, e47 (2019).
Article CAS PubMed PubMed Central Google Scholar
Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2010).
Article CAS PubMed Google Scholar
Robinson, M. D. & Oshlack, A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 11, R25 (2010).
Article PubMed PubMed Central Google Scholar
Robinson, M. D. & Smyth, G. K. Moderated statistical tests for assessing differences in tag abundance. Bioinformatics 23, 2881–2887 (2007).
Article CAS PubMed Google Scholar
McCarthy, D. J., Chen, Y. & Smyth, G. K. Differential expression analysis of multifactor RNA-seq experiments with respect to biological variation. Nucleic Acids Res. 40, 4288–4297 (2012).
Article CAS PubMed PubMed Central Google Scholar
Zhang, Y., Parmigiani, G. & Johnson, W. E. ComBat-seq: batch effect adjustment for RNA-seq count data. NAR Genom. Bioinform. 2, lqaa078 (2020).
Article PubMed PubMed Central Google Scholar
Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E. & Storey, J. D. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 28, 882–883 (2012).
Article CAS PubMed PubMed Central Google Scholar
Phipson, B., Lee, S., Majewski, I. J., Alexander, W. S. & Smyth, G. K. Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Ann. Appl Stat. 10, 946–963 (2016).
Article PubMed PubMed Central Google Scholar
Lun, A. T., Chen, Y. & Smyth, G. K. It’s DE-licious: a recipe for differential expression analyses of RNA-seq experiments using quasi-likelihood methods in edgeR. Methods Mol. Biol. 1418, 391–416 (2016).
Article PubMed Google Scholar
Lund, S. P., Nettleton, D., McCarthy, D. J. & Smyth, G. K. Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates. Stat. Appl. Genet. Mol. Biol. https://doi.org/10.1515/1544-6115.1826 (2012).
Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl Acad. Sci. USA 102, 15545–15550 (2005).
Article CAS PubMed PubMed Central Google Scholar
Luecken, M. D. & Theis, F. J. Current best practices in single-cell RNA-seq analysis: a tutorial. Mol. Syst. Biol. 15, e8746 (2019).
Article PubMed PubMed Central Google Scholar
Hafemeister, C. & Satija, R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 20, 296 (2019).
Article CAS PubMed PubMed Central Google Scholar
Becht, E. et al. Dimensionality reduction for visualizing single-cell data using UMAP. Nat. Biotechnol. 37, 38–44 (2019).
Article CAS Google Scholar
Finak, G. et al. MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 16, 278 (2015).
Article PubMed PubMed Central Google Scholar
Tran, H. T. N. et al. A benchmark of batch-effect correction methods for single-cell RNA sequencing data. Genome Biol. 21, 12 (2020).
Article CAS PubMed PubMed Central Google Scholar
Johnson, W. E., Li, C. & Rabinovic, A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8, 118–127 (2007).
Article PubMed Google Scholar

Download references

Acknowledgements

We thank S. Gregori and G. Amodio for help with neutrophil isolation and culture experiments; F. Di Salvo, F. Porzio and M. Tassara for patient recruitment and data management; the Center for Omics Sciences, the Flow cytometry Resource, Advanced Cytometry Technical Applications Laboratory, Centro Risorse Biologiche at Ospedale San Raffaele; and the Centro Universitario di Statistica per le Scienze Biomediche at Vita-Salute San Raffaele University. Figures were created with Adobe Illustrator and BioRender.com. V.B. and F.V.M. conducted this study as partial fulfillment of a PhD in Molecular Medicine (Basic and Applied Immunology and Oncology program) at Vita-Salute San Raffaele University. R.D.M. is a New York Stem Cell Foundation – Robertson Investigator. M.A.C. and P.S. are supported by grants from the Italian Association for Cancer Research (AIRC) (IG 20339) and the Italian Ministry of University and Research (PRIN 20177J4E75_004). A.A. is supported by the Italian Telethon Foundation (SR-Tiget grant award B02). E.M. and N.C. are supported by fellowships from Fondazione Umberto Veronesi. This study was supported by grants from the Italian Telethon Foundation (SR-Tiget grant award F04 to R.O.) and the Italian Ministry of Health (GR-201602362156 to R.O. and S.C.). Research in the R.O. laboratory is supported by the European Research Council (starting grant 759532, X-TAM) and by AIRC (MFAG 20247 and AIRC 5×1000 special program 22737).

Author information

These authors contributed equally. Elisa Montaldo, Eleonora Lusito, Valentina Bianchessi, Nicoletta Caronni.

Authors and Affiliations

San Raffaele Telethon Institute for Gene Therapy (SR-Tiget), IRCCS San Raffaele Scientific Institute, Milan, Italy
Elisa Montaldo, Eleonora Lusito, Valentina Bianchessi, Nicoletta Caronni, Serena Scala, Luca Basso-Ricci, Carla Cantaffa, Alice Masserdotti, Mattia Barilaro, Simona Barresi, Marco Genua, Francesco Maria Vittoria, Giulia Barbiera, Raffaella Di Micco, Andrea Ditadi, Matteo Maria Naldini, Bernhard Gentner, Luigi Naldini, Fabio Ciceri, Alessandro Aiuti & Renato Ostuni
Center for Omics Sciences (COSR), IRCCS San Raffaele Scientific Institute, Milan, Italy
Dejan Lazarevic
Hematology and Bone Marrow Transplantation Unit, IRCCS San Raffaele Scientific Institute, Milan, Italy
Carlo Messina, Elisabetta Xue, Sarah Marktel & Fabio Ciceri
Molecular Hematology Unit, IRCCS San Raffaele Scientific Institute, Milan, Italy
Cristina Tresoldi
Immunohematology and Transfusion Medicine Unit, IRCCS San Raffaele Scientific Institute, Milan, Italy
Raffaella Milani, Paola Ronchi, Salvatore Gattillo & Luca Santoleri
Pancreas Translational and Clinical Research Center, IRCCS San Raffaele Scientific Institute, Milan, Italy
Giulio Belfiori, Francesca Aleotti, Stefano Crippa & Massimo Falconi
Section of General Pathology, Department of Medicine, University of Verona, Verona, Italy
Elisa Gardiman, Nicola Tamassia, Marco Antonio Cassatella, Francesca Pettinella & Patrizia Scapini
Area of Cell and Developmental Biology, Centro Nacional de Investigaciones Cardiovasculares Carlos III, Madrid, Spain
Andrés Hidalgo
Singapore Immunology Network (SIgN), A*STAR, Biopolis, Singapore, Singapore
Immanuel Kwok & Lai Guan Ng
Vita-Salute San Raffaele University, Milan, Italy
Stefano Crippa, Massimo Falconi, Luigi Naldini, Fabio Ciceri, Alessandro Aiuti & Renato Ostuni
Pediatric Immunohematology and Bone Marrow Transplantation Unit, IRCCS San Raffaele Scientific Institute, Milan, Italy
Alessandro Aiuti

Authors

Elisa Montaldo
View author publications
You can also search for this author in PubMed Google Scholar
Eleonora Lusito
View author publications
You can also search for this author in PubMed Google Scholar
Valentina Bianchessi
View author publications
You can also search for this author in PubMed Google Scholar
Nicoletta Caronni
View author publications
You can also search for this author in PubMed Google Scholar
Serena Scala
View author publications
You can also search for this author in PubMed Google Scholar
Luca Basso-Ricci
View author publications
You can also search for this author in PubMed Google Scholar
Carla Cantaffa
View author publications
You can also search for this author in PubMed Google Scholar
Alice Masserdotti
View author publications
You can also search for this author in PubMed Google Scholar
Mattia Barilaro
View author publications
You can also search for this author in PubMed Google Scholar
Simona Barresi
View author publications
You can also search for this author in PubMed Google Scholar
Marco Genua
View author publications
You can also search for this author in PubMed Google Scholar
Francesco Maria Vittoria
View author publications
You can also search for this author in PubMed Google Scholar
Giulia Barbiera
View author publications
You can also search for this author in PubMed Google Scholar
Dejan Lazarevic
View author publications
You can also search for this author in PubMed Google Scholar
Carlo Messina
View author publications
You can also search for this author in PubMed Google Scholar
Elisabetta Xue
View author publications
You can also search for this author in PubMed Google Scholar
Sarah Marktel
View author publications
You can also search for this author in PubMed Google Scholar
Cristina Tresoldi
View author publications
You can also search for this author in PubMed Google Scholar
Raffaella Milani
View author publications
You can also search for this author in PubMed Google Scholar
Paola Ronchi
View author publications
You can also search for this author in PubMed Google Scholar
Salvatore Gattillo
View author publications
You can also search for this author in PubMed Google Scholar
Luca Santoleri
View author publications
You can also search for this author in PubMed Google Scholar
Raffaella Di Micco
View author publications
You can also search for this author in PubMed Google Scholar
Andrea Ditadi
View author publications
You can also search for this author in PubMed Google Scholar
Giulio Belfiori
View author publications
You can also search for this author in PubMed Google Scholar
Francesca Aleotti
View author publications
You can also search for this author in PubMed Google Scholar
Matteo Maria Naldini
View author publications
You can also search for this author in PubMed Google Scholar
Bernhard Gentner
View author publications
You can also search for this author in PubMed Google Scholar
Elisa Gardiman
View author publications
You can also search for this author in PubMed Google Scholar
Nicola Tamassia
View author publications
You can also search for this author in PubMed Google Scholar
Marco Antonio Cassatella
View author publications
You can also search for this author in PubMed Google Scholar
Andrés Hidalgo
View author publications
You can also search for this author in PubMed Google Scholar
Immanuel Kwok
View author publications
You can also search for this author in PubMed Google Scholar
Lai Guan Ng
View author publications
You can also search for this author in PubMed Google Scholar
Stefano Crippa
View author publications
You can also search for this author in PubMed Google Scholar
Massimo Falconi
View author publications
You can also search for this author in PubMed Google Scholar
Francesca Pettinella
View author publications
You can also search for this author in PubMed Google Scholar
Patrizia Scapini
View author publications
You can also search for this author in PubMed Google Scholar
Luigi Naldini
View author publications
You can also search for this author in PubMed Google Scholar
Fabio Ciceri
View author publications
You can also search for this author in PubMed Google Scholar
Alessandro Aiuti
View author publications
You can also search for this author in PubMed Google Scholar
Renato Ostuni
View author publications
You can also search for this author in PubMed Google Scholar

Contributions

E.M. and E.L. contributed to the design of the study, performed experiments, analyzed data, prepared figures and edited the manuscript; E.M. performed immunophenotype analyses together with S.S. and L.B.R. E.M., V.B. and N.C. generated transcriptome and functional data together with C.C., A.M., M.B., S.B., M.G., F.M.V., D.L., E.G., N.T. and F.P. E.L. performed computational and statistical analyses with help from G. Barbiera. C.M., E.X., S.M., C.T., R.M., P.R., S.G., L.S., G. Belfiori and F.A. selected and recruited study participants, collected biological samples and acquired clinical information under the supervision of S.C., M.F. and F.C. R.D.M., A.D., M.M.N., B.G., A.H., I.K., L.G.N. and L.N. provided key scientific inputs. M.A.C. and P.S. supervised E.G., N.T. and F.P. and contributed to data interpretation. A.A. supervised S.S. and L.B.R. and contributed to data interpretation. R.O. conceptualized and coordinated the study, acquired funding, analyzed the data and wrote the paper. All authors read and edited the manuscript.

Corresponding authors

Correspondence to Elisa Montaldo or Renato Ostuni.

Ethics declarations

Competing interests

The authors declare no competing interests

Peer review

Peer review information

Nature Immunology thanks Zvi Fridlender, Hongbo Luo and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available. Primary Handling Editor: L. A. Dempsey, in collaboration with the Nature Immunology team.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data

Extended Data Fig. 1 Leukocyte dynamics in G-CSF-treated donors.

a. Full gating strategy used to identify leukocyte subsets in whole PB or BM samples (blue: myeloid cells; green: lymphoid cells; red: HSPC; pink: CD45⁺ Lineage⁻CD34⁻cells; brown: CD45⁻ cells). Neutrophil subsets are numbered from 1 to 4 (1: SSC^hi CD38⁺ CD11c⁻ CD10⁻ neutrophil precursors; 2: SSC^hi CD38⁻ CD11c⁻ CD10⁻ immature neutrophils; 3: SSC^hi CD38⁻ CD11c⁺ CD10⁻ immature neutrophils; 4: SSC^hi CD38⁻ CD11c⁻ CD10⁺ mature neutrophils). b. Gating strategy used to identify low-density neutrophils (LDNs) within PBMCs. c. Gating strategy used to identify normal-density neutrophils (NDNs) within granulocytes. d. Representative May-Grunwald Giemsa staining of NDNs isolated by bead sorting (lower row) and LDN (upper row) isolated by FACS sorting from PB of healthy donors (steady state, n = 3) and from G-CSF mobilized PB (n = 2). e,f. Absolute HSPC counts in whole PB (e) and frequency of HSPC subsets (gated on Lin⁻ CD45⁺ CD34⁺ cells) in whole BM or PB (f) of controls or G-CSF-treated donors (PB n = 15; BM n = 14). g-j. Absolute counts of phenotypically defined hematopoietic stem cells (HSC), multipotent progenitors (MPP) and multi-lymphoid progenitors (MLP) (g), of committed common myeloid progenitors (CMP) or granulocyte-monocyte progenitors (GMP) (h) and of differentiated myeloid and lymphoid cells (i,j) in whole PB of controls or G-CSF-treated donors (n = 15). k. Representative contour plot showing immature and mature neutrophils in the whole PB of a representative G-CSF-treated donor. Gating strategies for the indicated cell types are also reported in Supplementary Table 7. Bar plot report data as mean ± SD. Numbers in red represent fold increases in the indicated conditions. Statistical analyses. e, g-j: two-sided Mann-Whitney test.

Source data

Extended Data Fig. 2 Leukocyte dynamics in HSC-T and PDAC patients.

a. Expression of the indicated markers in NDNs and LDNs from controls (n > 10) or G-CSF-treated donors (n > 12). b, c. Histogram (b) and cumulative histogram (c) plot showing the expression levels of CD49d in neutrophil precursors, immature and mature neutrophils. d, e. Contour (d) and cumulative histogram (e) plots showing percentages of EdU⁺ cells within NDNs and LDNs from BM, CB and PB samples of the indicated patients (BM n = 3, CB n = 6, G-CSF n = 8, HSC-T 2° f.u. n = 3). f, g. Representative histogram (f) and cumulative histogram (g) plot showing expression of the indicated markers in in NDNs and LDNs from controls (n > 11) or G-CSF-treated donors (n > 11). h. Representation of leukocytes dynamics in HSC-T patients. i. Quantification of white blood cells (WBCs), neutrophil, monocyte and lymphocyte count in HSC-T patients. Gray intervals highlight normal ranges. Arrows indicate the beginning of myeloablative conditioning; day 0 indicates the day of HSC-T. j. Counts of monocytes or lymphocytes in PB of controls (n = 8) or HSC-T patients (1° f.u. n = 8, 2° f.u. n = 9, 3° f.u. n = 3). k,l. Percentage of neutrophil precursors, immature and mature neutrophils within PBMCs (k) or LDNs (l) in controls (n = 8) and HSC-T patients (1° f.u. n = 7, 2° f.u. n = 8). m-o. Absolute counts (m) and frequencies of HSPC subsets (n, o) in whole PB (n = 15) or BM (n = 14) of controls and PB of PDAC patients (n = 8). p. Leukocyte counts in whole PB of controls (n = 15) and PDAC patients (n = 8). q. Neutrophil-to-lymphocyte ratio (NLR) in whole PB of controls (n = 15) and PDAC (n = 8) patients. NLR is calculated as the ratio between absolute counts (FACS) of neutrophils and total lymphocytes. r. Leukocyte counts (hemocytometer) and corresponding NLR values in whole PB of controls (n = 10), IPMN (n = 12) and PDAC (n = 15) patients. s. Contour plots showing CD16 and CD11b expression in LDNs of three PDAC patients. t, u. Percentage of neutrophil precursors, immature and mature neutrophils within LDNs (t) and PBMCs (u) of controls (n = 12), IPMN (n = 12) and PDAC (n = 18) patients. v. Percentage of EdU⁺ cells within neutrophil precursors, immature and mature neutrophils in PB of PDAC patients (n = 6). Gating strategies for the indicated cell types are reported in Supplementary Table 7. Bar plots report data as mean ± SD. Statistical analyses. a, c, g, j, k and r: Kruskal-Wallis test plus two-sided Dunn’s multiple comparison. m, o-q: two-sided Mann-Whitney test.

Source data

Extended Data Fig. 3 Purity of isolated cell populations.

a-d. Representative contour plots showing cell purity before and after magnetic bead selection of LDNs (a, b), monocytes (c) and NDNs (d) from PB samples. e. Representative contour plots showing the gating strategy used to isolate LDNs and monocytes (left panel) and post-sort purity analysis of sorted cells (right panel). f. Representative contour plots showing the gating strategy used to isolate BM neutrophil subsets and post-sort purity analysis of sorted cells.

Extended Data Fig. 4 Bulk RNA-Seq analysis of NDNs, LDNs and monocytes.

a, b. Principal-component analysis (PCA) plots of bulk RNA-Seq datasets of NDNs, LDNs and monocytes isolated from PB of healthy controls (n = 19), G-CSF-treated donors (n = 17), HSC-T (n = 8), PDAC (n = 15) and IPMN (n = 14) patients, as well as of neutrophil differentiation intermediates from BM of healthy donors (n = 3) and HSC-T patients (n = 7). Samples are colored based on cell type (a) or stress condition (b), as indicated by the legends. Filled area plot on the left show the frequency of neutrophil precursors (pre), immature (imm) and mature (mat) neutrophils for the corresponding NDNs and LDNs samples along PC2. c. PCA plots of bulk RNA-Seq datasets of NDNs, LDNs or monocytes. Colors represent stress condition, while shapes reflect the tissue of origin (PB circle; BM triangle), as indicated in the legend.

Extended Data Fig. 5 Validation of RNA-Seq analyses in NDNs.

a, b. Representative contour plots (a) and cumulative bar plot (b) showing the basal expression of IL-1b in NDNs and LDNs isolated from controls (n = 3) and G-CSF-treated donors (n = 3). c-e. Cumulative bar plots showing the expression of the indicated genes in NDNs isolated from controls and G-CSF-treated donors (c), HSC-T patients (d) or IPMN and PDAC patients (e). f,g. Representative histogram plots (f) and cumulative histogram plots (g) showing the expression of the indicated markers in NDNs isolated from PB of controls (n > 5) and G-CSF-treated donors (n > 5). Bar plots report data as mean ± SD. Statistical analyses. c and g: two-sided Mann-Whitney test. b, d and e Kruskal-Wallis test plus two-sided Dunn’s multiple comparison.

Source data

Extended Data Fig. 6 Plasma factors in G-CSF-treated donors, HSC-T or PDAC patients.

a,b. Concentration of selected factors in the plasma of controls (n = 19) and G-CSF-treated donors (n = 13) (a) or controls (n = 19) and IPMN (n = 15) or PDAC (n = 18) patients (b). c. Correlation between plasma concentrations of the indicated factor and frequencies of neutrophils precursors or LDNs in the PMBC fraction. Colors indicate calculated Spearman’s correlation coefficients (p-value < 0.05). Gray, not significant. Data are shown for all experimental conditions (upper heat map) or excluding G-CSF-treated donors (lower heat map) (steady state n = 14; G-CSF n = 9; HSC-T 1° f.u. n = 7; HSC-T 2° f.u. n = 8; IPMN n = 14; PDAC n = 16) d, e. Correlation between plasma concentrations of the indicated factor and frequencies of LDNs in the PMBC fraction combining all samples together (d) or excluding (e) G-CSF-treated donors. Spearman’s correlation and p-values are shown for each plot. Bar plots report data as mean ± SD Statistical analyses. a: two-sided Mann-Whitney test; b: Kruskal-Wallis test plus two-sided Dunn’s multiple comparison; c-e: Spearman’s correlation.

Source data

Extended Data Fig. 7 Single-cell RNA-Seq analyses of human neutrophils.

a. Gating strategy used to isolate CD15⁺ neutrophils from whole BM (upper panel) or PB (lower panel) samples. Expression of CD16 and CD11b in sorted cells is shown. b. UMAP plot showing donor or patient identities. c-d. UMAP plots showing the expression of gene modules related to neutrophil maturation identified in the indicated studies. e. UMAP plots showing the expression of gene modules identified from bulk RNA-Seq analysis (see Fig. 4a and Supplementary Table 9). f, g. Stacked bar plots showing the frequency of cells from PB or BM samples (f) or from donors and patients (g) for each neutrophil cluster. h. Model depicting divergent developmental trajectories in stress-elicited neutrophils, leading to diverse gene expression programs of mature cells.

Extended Data Fig. 8 CellHarmony analyses of neutrophils from G-CSF-treated donors, HSC-T or PDAC patients.

a. Heat map showing standardized average expression (computed on normalized expression levels) of developmental marker genes identified by CellHarmony and expressed in at least 20% of cells from reference datasets for the indicated neutrophil subsets in controls (reference, white bars) and PDAC patients (query, black bars). Color bars represent stages of neutrophil development after alignment of scRNA-Seq data with CellHarmony. The number of cells from reference and query datasets for each cluster is shown at the top; the number of developmental marker genes for each cluster is shown on the left. Selected representative genes are highlighted on the right. b, c. Filled area plots showing mean expression in scRNA-Seq data of selected developmental marker genes in neutrophil subsets from controls (gray) and PDAC (green) (b) or controls and G-CSF-treated donors (dark blue) or HSC-T patients (light blue) (c). Numbers on the x axis indicate the stages of neutrophil development identified by CellHarmony. d. Box plots showing standardized average expression of genes upregulated (see Methods) in the indicated neutrophil subsets from PDAC patients versus controls. Each plot refers to induced genes in query versus reference scRNA-Seq datasets for neutrophils at each stage of development defined by CellHarmony. Box plots represent the median, interquartile range (IQR), minimum (25th percentile, 1.5 × IQR) and maximum (75th percentile, 1.5 × IQR). Sample size corresponds to the number of cells indicated in the heat map (a). e. Venn diagram showing the overlap between upregulated genes in G-CSF-treated donors and HSC-T and PDAC patients in the indicated stages of neutrophil development. A selection of genes upregulated in all conditions and of stress-specific genes are shown (see Supplementary Table 37).

Extended Data Fig. 9 CellHarmony analyses of neutrophils from COVID-19 patients.

a. Scheme depicting the experimental and computational strategies used to isolate and process for scRNA-seq analysis cells from controls and COVID-19 patients. b. Heat map showing standardized average expression (computed on normalized expression levels) of developmental marker genes identified by CellHarmony and expressed in at least 20% of cells from reference datasets for the indicated neutrophil subsets in controls (reference, white bars) and COVID-19 patients (query, black bars). Color bars represent stages of neutrophil development after alignment of scRNA-Seq data with CellHarmony. The number of cells from reference and query datasets for each cluster is shown at the top; the number of developmental marker genes for each cluster is shown on the left. Selected representative genes are highlighted on the right. c. Box plots showing standardized average expression of genes up regulated (see Methods) in the indicated neutrophil subsets from COVID-19 patients versus controls. Each plot refers to induced genes in query versus reference scRNA-Seq datasets for neutrophils at each stage of development defined by CellHarmony. Box plots represent the median, interquartile range (IQR), minimum (25th percentile, 1.5 × IQR) and maximum (75th percentile, 1.5 × IQR). Sample size corresponds to the number of cells indicated in the heat map (b). d. Bar plots showing NES of selected GO categories enriched within genes expressed at higher levels in neutrophil subsets from COVID-19 patients as compared to controls. Colors represent stages of neutrophil development defined by CellHarmony. e. Violin plots showing normalized expression levels of selected genes induced in mature neutrophils from COVID-19 patients as compared to controls. Colors represent stages of neutrophil development defined by CellHarmony.

Extended Data Fig. 10 Single-cell RNA-Seq analysis of IFN-stimulated neutrophils.

a. Experimental strategy used to enrich and mix LDNs and NDNs from CB samples, ensuring a sufficient representation of all neutrophil subsets (see Methods). b, c. Representative contour plots (b) and stacked bar plot (c) showing the percentage of neutrophil precursors (pre), immature (imm) and mature (mat) neutrophils in LDN, NDN and LDN-NDN mix (1:3). d. Schematic representation of CB neutrophil stimulation and processing for scRNA-Seq analysis. e. tSNE plots showing the expression of selected developmental marker genes. f. tSNE plots showing expression of cluster 3 marker genes (corresponding to CB neutrophils) in PB and BM neutrophils from steady-state controls.

Source data

Supplementary information

Reporting Summary

Peer Review File

Supplementary Table 1

Supplementary Tables 1–44.

Source data

Source Data Fig. 1

Statistical source data.

Source Data Fig. 2

Statistical source data.

Source Data Fig. 3

Statistical source data.

Source Data Fig. 5

Statistical source data.

Source Data Fig. 6

Statistical source data.

Source Data Extended Data Fig. 1

Unprocessed microscope images.

Source Data Extended Data Fig. 1

Statistical source data.

Source Data Extended Data Fig. 2

Statistical source data.

Source Data Extended Data Fig. 5

Statistical source data.

Source Data Extended Data Fig. 6

Statistical source data.

Source Data Extended Data Fig. 10

Statistical source data.

Rights and permissions

Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions

About this article

Cite this article

Montaldo, E., Lusito, E., Bianchessi, V. et al. Cellular and transcriptional dynamics of human neutrophils at steady state and upon stress. Nat Immunol 23, 1470–1483 (2022). https://doi.org/10.1038/s41590-022-01311-1

Download citation

Received: 26 November 2021
Accepted: 10 August 2022
Published: 22 September 2022
Issue Date: October 2022
DOI: https://doi.org/10.1038/s41590-022-01311-1

This article is cited by

Quantitative proteomics reveals tissue-specific, infection-induced and species-specific neutrophil protein signatures
- Gabriel Sollberger
- Alejandro J. Brenes
- Andrew J. M. Howden
Scientific Reports (2024)
Neutrophils: from IBD to the gut microbiota
- Camille Danne
- Jurate Skerniskyte
- Harry Sokol
Nature Reviews Gastroenterology & Hepatology (2024)
Patients with schizophrenia and bipolar disorder display a similar global gene expression signature in whole blood that reflects elevated proportion of immature neutrophil cells with association to lipid changes
- Anja Torsvik
- Hans-Richard Brattbakk
- Vidar M. Steen
Translational Psychiatry (2023)
Charting granulopoietic disturbances in sepsis
- Roza Maria Barouni
- Renato Ostuni
Nature Immunology (2023)
Strategies of neutrophil diversification
- Miguel Palomino-Segura
- Jon Sicilia
- Andrés Hidalgo
Nature Immunology (2023)

Subjects

Abstract

Similar content being viewed by others

Main

Results

Dynamics and phenotype of G-CSF-elicited human neutrophils

Dynamics and phenotype of neutrophils during HSC-T or PDAC

Functional properties of G-CSF-elicited neutrophils

Bulk transcriptome analysis of monocytes, NDNs and LDNs

Transcriptional responses to stress in NDNs and monocytes

Analysis of plasma factors elicited by G-CSF, HSC-T or PDAC

Transcriptional diversity of human neutrophils upon stress

Dynamics of neutrophil differentiation during stress

Transcriptional changes of human neutrophils during stress

Transcriptional responses to IFN by developing neutrophils

Discussion

Methods

Experimental methods

Study participants and sample collection

Controls and G-CSF-treated donors

Patients receiving HSC-T

Patients with PDAC and IPMN

Cell isolation

Flow cytometry

Whole blood staining

Mononuclear cells and granulocyte staining

Fluorescence activated cell sorting

EdU incorporation

ROS production

NETosis

Ex vivo stimulation of NDNs and LDNs

Plasma collection

ELISA

Cytospin and May-Grünwald Giemsa staining

Real-time quantitative PCR

Ex vivo stimulation of CB neutrophils

Bulk RNA sequencing

Single-cell RNA sequencing

Computational methods

Bulk RNA-seq analyses on NDNs, LDNs and monocytes

Data processing

Heat map of variable genes

Differential gene expression analysis

Principal-component analysis

Gene set enrichment analyses

Single-cell RNA-seq analyses of PB and BM neutrophils

Data processing

Gene set enrichment analysis

cellHarmony analyses

Gene set enrichment analysis

Single-cell RNA-seq analyses of CB neutrophils

Statistics and reproducibility

Reporting summary

Data availability

References

Acknowledgements

Author information

Authors and Affiliations

Contributions

Corresponding authors

Ethics declarations

Competing interests

Peer review

Peer review information

Additional information

Extended data

Supplementary information

Source data

Rights and permissions

About this article

Cite this article

Share this article

This article is cited by

Search

Quick links