Experimental methodsStudy 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−1 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−1 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 NH4Cl, 10 mM KHCO3 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 described19. 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 previously17. 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−1. 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 described17 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+CD16int and BM4 CD11b+CD16hiCD10+ 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 106 cells ml−1 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 × 106 cells ml−1 in Hank’s Balanced Salt Solution (HBSS), pH 7.4, supplemented with 10% FBS, 0.5 mM CaCl2 and 1 mg ml−1 glucose. O2− production in response to 20 ng ml−1 PMA (Sigma) was assessed by the Cytochrome C reduction assay (Cayman), as previously described51. For flow cytometry analysis of ROS, 1 × 105 mononuclear cells or granulocytes were incubated with dihydrorhodamine-123 (Cayman Chemical) and left untreated or stimulated with PMA 20 ng ml−1 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 × 106 cells ml−1 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 × 106 ml−1 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 106 cells ml−1 in RPMI 1640 + 10% FBS + 1% Gln + 1% pen/strep alone or with G-CSF, IFN-β or IFN-γ all used at 10 ng ml−1. 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 method52 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 methodsBulk 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)53. The FeatureCounts function from Rsubread package (v.3.7)54 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)55. 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 method56 using calcNormFactors function. The weighted likelihood empirical Bayes method57,58 was used to calculate the posterior dispersion estimates through the estimateDisp function. The ComBat_seq function59 from the sva package (v.3.38.0)60 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

Log2 (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 estimation57,61 using the glmQLFit function and after computing the dispersion with estimateDisp function. A quasi-likelihood F-test62,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(log2FC) ≥ 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 log2FC in stress versus steady state. Gene set enrichment analyses (GSEA) (v.4.0.3)64 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 count65. We also applied a gene-wise filter on the average counts to remove low-abundance genes62. 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 normalization66 (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 UMAP67. Cluster-specific marker genes were identified using the MAST method68 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-Schrepping27 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 log2FC > 0. GSEA (v.4.0.3)64 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 ComBat69,70 from the sva package (v.3.38.0) to adjust for potential batch effects between donors of the same condition. cellHarmony26 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 log2FC ≥ 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 log2FC ≥ 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 log2FC 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.