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

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.

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

a. Full gating strategy used to identify leukocytes 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.

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, HSCT 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.

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 HSCT 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-1β 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.

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 heatmap) or excluding G-CSF-treated donors (lower heatmap) (steady state n= 14; G-CSF n=9; HSCT 1° f.u. n=7; HSCT 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.

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. Heatmap showing standardized average expression (computed on normalized expression levels) of developmental marker genes identified by Cell Harmony 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 Cell Harmony. 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 (grey) 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 Cell Harmony. d. Box plots showing standardized average expression of genes up regulated (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 Cell Harmony. 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 heatmap (a). e. Venn diagram showing the overlap between up-regulated genes in G-CSF treated donors and HSC-T and PDAC patients in the indicated stages of neutrophil development. A selection of genes up-regulated 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 scRNAseq analysis cells from controls and COVID-19 patients. b. Heatmap showing standardized average expression (computed on normalized expression levels) of developmental marker genes identified by Cell Harmony 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 Cell Harmony. 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 Cell Harmony. 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 heatmap (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 Cell Harmony. 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 Cell Harmony.

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.

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

a. Absolute counts of neutrophils (gated as CD45⁺ CD33⁺/CD66b⁺ SSC^hi) in whole PB of controls or G-CSF-treated donors (n=15). b, c Absolute counts (b) and percentage (c) of neutrophil subsets (gated as shown in Extended Data Fig. 1a) in whole PB of controls (n=15) or G-CSF-treated donors (n=15). d, e. Representative contour plots (d) or cumulative histogram plots (e) showing the frequencies of low-density neutrophils (LDNs) in the PBMC fraction of controls (n=16) and G-CSF-treated donors (n=15). f. Histogram plots showing representative expression levels of the indicated markers in normal density neutrophils (NDNs) from PB of controls (n>10) or G-CSF-treated donors and LDN of G-CSF-treated donors (n>12). g. Gating strategy used to identify neutrophil precursors (pre), immature (imm) and mature (mat) neutrophils within LDNs. h, i. Percentage of neutrophil precursors, immature and mature neutrophils within the PBMC fraction (h) or within LDNs (i) in controls (n=12) and G-CSF-treated donors (n=17). j, k. Representative contour plots (j) and cumulative histogram plots (k) showing percentages of EdU⁺ cells within neutrophil precursors, immature and mature neutrophils in controls (n=3) or G-CSF-treated donors (n=8). Gating strategies for the indicated cell types are reported in Extended Data Fig. 1a-c and Supplementary Table 7. Bar plots represent data as mean ± SD. Numbers in red represent fold increases in the indicated conditions. Statistical analyses. a, b, e, and h: two-sided Mann-Whitney test.

Fig. 2. Dynamics and phenotype of neutrophils during HSC-T or PDAC.

a. Absolute counts of neutrophils in whole PB of controls (n=8) or HSC-T patients (1° f.u. n=8, 2° f.u. n=9, 3° f.u. n=3). b, c. Absolute counts (b) and percentage (c) of neutrophil subsets (gated as shown in Extended Data Fig. 1a) in whole PB of controls (n=8) or HSC-T patients (1° f.u. n=7, 2° f.u. n=8, 3° f.u. n=3). d, e. Representative contour plots (d) or cumulative histogram plots (e) showing the frequencies of LDNs in the PBMC fraction of controls (n=16) and HSC-T patients (1° f.u. n=8, 2° f.u. n=8, 3° f.u. n=5). f, g. Representative contour plots (f) and cumulative histogram plots (g) showing percentages of EdU⁺ cells within neutrophil precursors, immature and mature neutrophils in PB samples of HSC-T patients (2° f.u. n=3). h, i. Representative contour plots (h) or cumulative histogram plots (i) showing the frequencies of LDNs in the PBMC fraction of controls (n=16), IPMN (n=14), and PDAC patients (n=16). j. Representative contour plots showing percentages of EdU⁺ cells within neutrophil precursors, immature and mature neutrophils in PB samples of PDAC patients. Gating strategies for the indicated cell types are reported in Extended Data Fig. 1a-c and Supplementary Table 7. Bar plots report data as mean ± SD. Statistical analyses. a, b, e and i: Kruskal-Wallis test plus two-sided Dunn’s multiple comparison.

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

a. Schematic description of ex vivo experiments performed to evaluate ROS production, NETosis, and cytokine release by neutrophils. b. Representative histogram plots showing Rhodamine 123 signal in PMA-stimulated NDNs and LDNs from controls and G-CSF treated donors. c. Line plot showing percentage of ROS⁺ cells in PMA-stimulated NDNs and LDNs from controls (n=2) and G-CSF treated donors (n=5). d. Representative histogram plots showing Rhodamine 123 signal in PMA-stimulated neutrophil precursors, immature, and mature neutrophils from G-CSF-treated donors. e. Line plot showing percentage of ROS⁺ cells in PMA-stimulated neutrophil precursors, immature, and mature neutrophils from G-CSF-treated donors (n=5). f. Line plots showing ROS levels in PMA-stimulated NDNs and LDNs from controls (n=3) or G-CSF-treated donors (n=3). g. Cumulative histogram plot showing PMA-induced NET release in NDNs and LDNs from controls (n=5) and G-CSF treated donors (n=6). h. Cumulative histogram plots showing the concentration of indicated cytokines released by R848-stimulated NDNs and LDNs from controls (n=6) or G-CSF-treated donors (n=8). Gating strategies for the indicated cell types are reported in Supplementary Table 7. Bar plots and line charts report data as mean ± SEM. Statistical analyses. c, e, f, g and h: Two-way ANOVA plus Tukey's multiple comparisons test. In panel f asterisks refer to the comparison between NDNs from controls versus G-CSF-treated donors; hashes refer to the comparison NDNs from controls versus LDNs from G-CSF-treated donors. * or # p-value < 0.05; ** or ## p-value < 0.01; *** or ### p-value < 0.001; **** or #### p-value < 0.0001; full p-values are reported in source data.

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

a. Heatmap showing normalized expression levels (Z-score) of variable genes (n=1,684, see Methods) in NDNs, LDNs and monocytes isolated from PB or BM 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 developing BM neutrophils from healthy donors (n=3) and HSCT patients (n=7). The row dendrogram represents hierarchical clustering of gene modules identified by k-means, and the column dendrogram represents hierarchical clustering of RNA-Seq samples. Legends and color bars at the top indicate sample identities by cell type and at the bottom by experimental condition. Numbers on the row dendrogram represent the identity and size of each gene module, with representative transcripts shown on the right. Gating strategies for cell sorting are reported in Extended Data Fig.3a-f and Supplementary Table 7. See Supplementary Table 8 for the full list of samples (n=210). b. Volcano plots showing differentially expressed genes in NDNs from G-CSF-treated donors, HSC-T (1^st and 2^nd follow-up) or PDAC patients as compared to steady-state controls. The x- and y-axes indicate the expression fold change (log2) and the false discovery rate (FDR) (-log₁₀) for each gene versus controls, respectively. Legends highlight up-regulated (red) or down-regulated (blue) transcripts, as well as genes not passing cut-off criteria for fold change (black) and FDR (grey) (see Methods). Selected representative genes are shown. c. Venn diagram showing the overlap between genes up-regulated in NDNs isolated form G-CSF treated donors or from HSC-T patients. Genes in green are also induced in NDNs isolated from PDAC patients. See Supplementary Table 13. d. Bar plots showing the normalized enrichment score (NES) of selected gene ontology (GO) categories enriched within genes up-regulated (red) or down-regulated (blue) in NDNs from the indicated experimental condition versus controls. e. Bar plot showing the number of genes induced (log₂FC >1.5 and FDR < 0.05) in NDNs and monocytes isolated from G-CSF-treated donors or from HSC-T patients.

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

a. Bar plots showing the fold change of the mean concentration of the indicated factors in the plasma of G-CSF-treated donors, HSC-T, IPMN or PDAC patients as compared to controls (steady state n= 19; G-CSF n=13; HSCT 1° f.u. n=9; HSCT 2° f.u. n=9; IPMN n=15; PDAC n=18). b, c. Concentration of selected factors in the plasma of controls and G-CSF-treated donors (b) or controls and HSC-T patients at the indicated follow-ups (c) (steady state n=19; G-CSF n=13; HSCT 1° f.u. n=9; HSCT 2° f.u. n=9). d, e. Correlation between plasma concentrations of the indicated factors and frequencies of neutrophil precursors 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 (steady state n= 14; G-CSF n=9; HSCT 1° f.u. n=7; HSCT 2° f.u. n=8; IPMN n=14; PDAC n=16). Cumulative bar plots report data as mean ± SD. Statistical analyses. a: Wilcoxon signed-rank test followed by FDR calculation with two-stage step-up method of Benjamini, Krieger and Yekutieli; b: two-sided Mann-Whitney test; c: Kruskal-Wallis test plus two-sided Dunn’s multiple comparison; d and e: Spearman’s correlation. a: *p-value < 0.05; **p-value < 0.01; *** p-value < 0.001; ****p-value < 0.0001; full p-values are reported in source data.

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

a. UMAP plot showing scRNA-Seq transcriptomes of 130,628 cells, sorted as CD15⁺ neutrophils from whole PB or BM samples of healthy controls (PB n=2, BM n=2), G-CSF-treated donors (n=4), HSC-T (PB n=3, BM n=2) and PDAC (n=5) patients. Colors and numbers indicate clusters at resolution 1.5. Representative marker genes are shown for selected clusters. Groups of clusters corresponding to developing neutrophil subsets are indicated on the right. Pie charts report the frequency of PB or BM cells and of cells from controls, G-CSF-treated donors, HSC-T and PDAC patients. Clusters 25, 26 and 27 were classified as contaminants. b. UMAP plots showing colored according to tissue of origin (PB/BM) and the stress condition. c. Heatmap showing expression of up to 50 marker genes for each scRNA-Seq cluster, with selected transcripts highlighted on the left or on the right. Color bars indicate cluster identities. Clusters of contaminants are not shown. d. Bar plots showing NES of selected GO categories enriched within combined marker genes of clusters 4 and 7 (corresponding to mature PB neutrophils from HSC-T patients at the 1^st follow-up). e, f. Box plots showing the expression levels of combined marker genes of clusters 4 and 7 in neutrophil precursors, early immature, immature, and mature neutrophils from PB (e) or BM (f) samples from steady-state controls and HSC-T patients at the indicated follow-ups. Sample sizes are reported in Supplementary Table 27. FDR-adjusted p-values were calculated by two-sided Wilcoxon rank-sum test. Box plots represent the median, interquartile range (IQR), minimum (25th percentile, 1.5 × IQR) and maximum (75th percentile, 1.5 × IQR).

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

a, b. Heatmaps showing standardized average expression (computed on normalized expression levels) of developmental marker genes identified by Cell Harmony and expressed in at least 20% of cells from reference datasets, for the indicated neutrophil subsets in controls (reference, white bars) and G-CSF-treated donors (a) or HSC-T patients (b) (query, black bars). The following samples were included in the Cell Harmony analysis: PB and BM for healthy controls; PB (all time points) and BM (day 30 and > 180 post-transplant) for HSC-T patients; PB from G-CSF-treated donors; PB from PDAC patients. Colored bars and numbers represent stages of neutrophil development (1-precursors, 2-proliferating, 3-early immature, 4-immature, 5-mature BM, and 6-mature PB) after alignment of scRNA-Seq data with Cell Harmony (see Methods). 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. Filled area plots showing mean expression in scRNA-Seq data of selected developmental marker genes in neutrophil subsets from controls (grey), G-CSF-treated donors (dark blue) or HSC-T patients (light blue). Numbers on the x-axis indicate the stages of neutrophil development identified by Cell Harmony. d, e. Box plots showing standardized average expression of genes up-regulated (see Methods) in the indicated neutrophil subsets from HSC-T patients (d) or G-CSF-treated donors (e) versus controls. Each plot refers to induced genes in query versus reference scRNA-Seq datasets for neutrophils at each stage of development defined by Cell Harmony. 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 heatmaps (a, b). f, g. Violin plots showing normalized expression levels of selected genes induced in mature neutrophils from HSC-T patients (f) or G-CSF-treated donors (g) as compared to controls. Colors represent stages of neutrophil development defined by Cell Harmony. h, i. Bar plots showing NES of selected GO categories enriched within genes expressed at higher levels in neutrophil subsets from HSC-T patients (h) or G-CSF-treated donors (i) as compared to controls. Colors represent stages of neutrophil development defined by Cell Harmony.

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

a, b. t-distributed stochastic neighbor embedding (tSNE) plots showing scRNA-Seq transcriptomes of 22,240 cells, isolated as neutrophils (see Methods and Extended Data Fig. 10a-d) from CB samples and stimulated ex vivo with IFN-β or IFN-γ for 4 hours. Colors and numbers represent clusters at resolution 0.3 (a), or the type of treatment (b). Representative marker genes are shown for selected clusters. Groups of clusters corresponding to developing neutrophil subsets (precursors, immature and mature) are indicated. Clusters 7 and 8 were classified as contaminants. c. Stacked bar plots showing the fraction of cells corresponding to control (blue), IFN-β (red) or IFN-γ (green) treatment conditions for the indicated scRNA-Seq clusters. d. Stacked bar plots showing the fraction of cells belonging to scRNA-Seq clusters for the indicated culture condition. e. Heatmap showing normalized expression (Z-score) of up to 50 marker genes for each scRNA-Seq cluster, with selected transcripts highlighted on the right. Cluster identities and corresponding classifications as precursors, immature or mature neutrophils is shown by color bars at the bottom. Color bars at the top indicate cells corresponding to control (blue), IFN-β (red) or IFN-γ (green) treatments. f, g. Violin plots showing mean standardized expression of top 25 marker genes of cluster 4 (mature neutrophils, IFN-β-treated) (f) or cluster 5 (mature neutrophils, IFN-γ-treated) (g) in cells corresponding to neutrophil precursors, immature and mature neutrophils from controls or the indicated stimulation conditions; p-values were calculated by two-sided Wilcoxon rank-sum test. h, i. Violin plots showing normalized expression of selected genes induced by IFN-β (h) or IFN-γ (i) in cells corresponding to neutrophil precursors, immature and mature neutrophils from controls or the indicated stimulation conditions. j, k. Violin plots showing normalized expression of the genes encoding for IFN receptors and signaling molecule STAT1 in cells corresponding to neutrophil precursors, immature and mature neutrophils from controls or the indicated stimulation conditions.