Architecture of the neutrophil compartment

Mice

All experiments were performed on 6-to-24-week-old C57BL/6 male and female mice. Young mice were defined as 8 to 12 weeks old, and old mice were defined as 22 to 24 months old at the time of analysis. Mice were maintained under specific pathogen-free conditions with chow and water provided ad libitum. mouse lines used were on the C57B1/6 J background and housed under specific pathogen-free conditions at the Centro Nacional de Investigaciones Cardiovasculares Carlos III, Singapore Immunology Network or Yale University. All mouse husbandry and experimentation was conducted using protocols approved by local animal ethics committees and authorities. Mice (Mus musculus) were maintained in racks with individual ventilation cages according to current Spanish, Singapore and US legislation (RD 53/2013 and EU Directive 63/2010, respectively). Mice have access to dust- and pathogen-free bedding, as well as sufficient nesting and environmental enrichment materials, to facilitate nesting. All mice were kept in environmental conditions of 45–65% of relative humidity, temperature of 21–24 °C, and a light:dark cycle of 12 h:12 h. Mice with neutrophil-specific deficiency in Tgfbr2 (Tgfbr^ΔN) were generated by crossing MRP8^CRE mice with Tgfbr2^fl/fl mice⁴¹. Similarly, we generated neutrophil-specific mutants by crossing Junb-floxed⁴², Csf2r-floxed⁴³ and Ifnar1-floxed mice⁴⁴ with the MRP8^CRE driver. Apoe^–/– mice (B6.129P2-Apoetm1Unc; Taconic M&B). Ly6g^creERT2 mice were crossed with Rosa26^Tdtomato mice as in ref. ¹, resulting in the iLy6g^tdTom mice used in our fate mapping experiments. Gavage administration of tamoxifen (2 mg per mouse) was performed to induce CRE recombinase activity in 6-to-12-week-old male iLy6g^tdTom mice. JAXBoy (Ptprc^K302E) from Jackson laboratories and Tet2^−/− mice¹⁹ were used for adoptive bone marrow cell transfer. Eight-week-old male Germ-free mice (C57Bl/6) were kindly provided by the laboratory of N. Palm. In brief, Germ-free C57BL/6 mice were bred and maintained under sterile conditions in flexible film isolators (Class Biologically Clean) in the Palm laboratory Gnotobiotic Facility at Yale School of Medicine. Mice were housed in a temperature- and humidity-controlled room under a dark:light cycle of 12 h:12 h. All animal protocols were approved by the Yale University Institutional Animal Care and Use Committee (protocol 11513).

For the rewilding experiments, litters of mice were generated from multiple breeding pairs and randomly assigned to either remain in the institutional vivarium (laboratory mice) or be released into the outdoor enclosures (rewilded mice) to control for the microbiota at the onset of the experiment. Outdoor enclosures were previously described⁴⁵ and the protocols for releasing laboratory mice into the outdoor enclosure facility and then returning them to vivaria were approved by Princeton University (protocol 1982) and Rutgers University (protocol PROTO999900794). All protocols were approved by the corresponding local authorities of Madrid, Singapore, Rutgers, Princeton and Yale University.

Mouse models of disease

Stroke

Thrombotic occlusion of the middle cerebral artery was induced by the ferric chloride (FeCl₃) stroke model. In brief, mice were anaesthetized and maintained at 2% sevoflurane in a mixture of 0.2 l min⁻¹ O₂:0.8 l min⁻¹ air and temperature was kept at 36.5–37 °C using a heating blanket. The scalp was opened, and the middle cerebral artery was visualized with a stereomicroscope (PZMIV, World Precision Instruments). A piece of Whatman filter paper strip soaked in freshly prepared FeCl₃ (20%) was placed over the intact dura mater on the artery for 10 min and then removed to allow the formation of a thrombus. Following surgery, individual mice were returned to their cages with free access to water and food. Brains were collected 24 h after surgery to perform transcriptomics analysis.

Liver cholestasis

For the liver injury model, mice were fed a 0.1% of 3,5-diethoxycarbonyl-1,4-dihydrocollidine (DDC)-supplemented diet for 3 weeks before sample collection, housed with a 12 h:12 h light:dark cycle, and permitted ad libitum consumption of water as described⁴⁶. An additional group of mice was fed a 0.1% DDC-supplemented diet for three weeks and afterward allowed to recover for three days under standard mouse diet to study the reversibility of the cholestatic phenotype.

Influenza A infection

A stock of the virus strain A/PR8/34 (H1N1) was diluted, and 100 plaque-forming units were administered intranasally to isoflurane-anaesthetized 8-to-12-week-old male mice in 50 μl of PBS. Mouse weight was monitored daily after infection and mice that presented weight loss of more than 20% of their initial body weight were euthanized and considered deceased. For transcriptomic studies, blood and lungs were collected on day 4 after infection.

Pancreatitis

Acute pancreatitis was induced by intraperitoneal injections of 50 µg kg⁻¹ of cerulein (Sigma-Aldrich), every hour, for a total of 7 administrations. Mice were euthanized 24 h after the first injection.

Orthotopic pancreatic tumour model

Mice were anaesthetized with ketamine/xylazine, and had their abdomen shaved and swabbed with antiseptic. A 5 mm vertical incision was made in the skin and abdominal layer at a point 1 cm down from the xiphoid process of the sternum, and 1 cm to the right of the midline. The pancreas was exposed, 10⁵ FC1242 cells were resuspended in phosphate-buffered saline (PBS) and mixed with Matrigel (BD) in a 1:1 ratio and were injected as a volume of 50 µl into the body of the pancreas to form a visible bolus using a 30G insulin needle. The pancreas was then returned to the abdominal cavity. The abdominal layer was closed with absorbable 5/0 sutures, while the skin was closed with non-absorbable 5/0 sutures. Superglue was applied over the sutures to ensure that they did not come undone after surgery. Mice were resuscitated with saline and were subcutaneously administered Buprenorphine (10 mg kg⁻¹) and Enrofloxacin (Baytril, 1.5 mg kg⁻¹) for the 2 days following surgery. Mice were euthanized at week 5 following surgery and tissues were collected for transcriptomic analysis.

Orthotopic breast tumour model

Mice were implanted orthotopically with 5 × 10⁵ E0771 breast cancer cells in 50 μl Matrigel into the thoracic mammary gland of C57BL6/J mice. Additionally, the same procedure was followed using the BALB/c-derived 4T1 breast tumour cell cancer in BALB/c mice. Mice were euthanized at week 4 after implantation and tumours were collected for transcriptomic analysis.

Orthotopic lung cancer model

We administered 2 × 10⁵ LLC cells in 100 µl PBS intravenously into the lateral tail vein of 8-week-old C57BL6/J mice. Mice were euthanized at week 1, 2 or 3 after the implantation and bloods and lungs were used for transcriptomics analysis.

Subcutaneous lung cancer model

Mouse LLC1 implants were generated in 8-week-old C57BL/6 mice by subcutaneous implantation of 0.5 × 10⁶ cells (1 injection site per mouse). Tumour growth was followed every 2 days by measuring the 2 orthogonal external diameters using a calliper. Tumour volume was calculated as V = π/6 × L × W × H, where L, W and H represent length, width, and height, respectively. Tumours were excised and processed for flow cytometry analysis when they reached 0.5 cm³.

Caecal ligation and puncture-induced sepsis

Caecal ligation and puncture were performed as described⁴⁷. In brief, the peritoneal cavity of ketamine/xylazine-anaesthetized mice was exposed with a small incision and the caecum was exteriorized. 80% of the caecum distal to the ileo-caecal valve was ligated using non-absorbable 7-0 suture. A 23G needle was then used to puncture both walls of the distal end of the caecum, and a small drop of faeces was extruded through the perforation. The ligated and punctured caecum was relocated inside the peritoneal cavity and both peritoneum and skin were closed. Blood was extracted three days after the puncture.

Peritonitis

Male 8-to-12-week-old mice were injected intraperitoneally with zymosan (1 mg, intraperitoneal injection, Sigma). After 2 h and 72 h we performed a peritoneal lavage for transcriptomic studies.

Myocardial infarction

Male 8-to-12-week-old mice were intubated, and temperature controlled throughout the experiment at 36.5 °C to prevent hypothermic cardioprotection. Thoracotomy was then performed, and the left anterior descending artery was ligated with a nylon 8/0 monofilament suture for 45 min. At the end of the ischaemia, the chest was closed, and mice were kept with 100% O₂ and given analgesia with buprenorphine (subcutaneous injection, 0.1 mg kg⁻¹) as described previously⁴⁷. Mice were euthanized 24 h or 72 h h after surgery and the heart was isolated for transcriptomics studies.

Bleomycin-induced fibrosis

We administered 1 mg kg⁻¹ of bleomycin sulfate to 8-to-12-week-old mice as previously described⁴⁸. In brief, bleomycin was dissolved in saline and was instilled into the tracheal lumen through a cannula under isoflurane (2.5% in oxygen) anaesthesia. Bleomycin was injected at day 0 and at day 4. Mice were euthanized three weeks after bleomycin injection and lungs were collected for transcriptomics analysis.

Staphylococcus aureus Infection

Mice were intravenously infected with 2.5 × 10⁷ CFU of S. aureus (RNU4220 strain) and monitored for weight loss. For single-cell transcriptomic studies, blood was collected five days after infection.

Candida infection

Mice were intravenously infected with 1.5 × 10⁵Candida albicans conidia (SC5314 strain), blood for single-cell transcriptomic studies was extracted at day 6 after infection⁴⁷.

LPS-induced inflammation

For transcriptomic studies, 400 ng of LPS (Sigma) were injected intravenously. Blood and tissues were collected 24 h after injection. An intraperitoneal lethal dose of LPS (40 mg kg⁻¹) was used as a model of endotoxic shock. Mice were monitored daily for weight loss. A weight loss of more than 20% of initial body weight was considered a lethal event and mice were euthanized at humane endpoints.

High-fat diet

Apoe^−/− mice were fed for 6 weeks with a control or high-fat diet (HFD, 10.7% total fat, 0.75% cholesterol; Sniff) before blood extraction.

Dextran sulfate sodium colitis

To induce colitis, mice received for 9 days water with 1.5% dextran sulfate sodium salt (MP Biomedicals) as previously described⁴⁹. Blood was collected on day 9 after dextran sulfate sodium treatment.

Hindlimb ischaemia

Hindlimb ischaemia experiments were performed as described⁵⁰. In brief, mice were anaesthetized with isoflurane, the hindlimb was shaved, and, following a small incision in the skin, both the proximal end of the femoral artery and the distal portion of the saphenous artery were ligated. The artery and all side-branches were dissected free; after this, the femoral artery and attached side-branches were excised. Immediately after surgery, perfusion was measured by laser Doppler imaging of plantar regions of interest (ROIs) (Moor Instruments) and calculated as ratio of left (ligated) versus right (unligated) values. Ischaemic muscle samples for transcriptomics analysis were collected one day after surgery.

Model of clonal haematopoiesis and PCSK9-induced hypercholesterolemia

To model TET2 loss-of-function-driven clonal haematopoiesis, we performed an adoptive bone marrow transfer without pre-conditioning. Ten-week-old unirradiated JAXBoy (Ptprc^K302E) recipient mice were intravenously injected with a total of 1.5 × 10⁷ unfractionated CD45.2⁺ Tet2^−/− bone marrow cells, administered as 3 consecutive daily doses of 5 × 10⁶ cells⁵¹. Donor cells were collected from age-matched littermate Tet2^−/− mice (8 to 10 weeks old) by flushing femurs and tibias following euthanasia. To induce hypercholesterolaemia, a recombinant AAV vector encoding a gain-of-function form of mouse PCSK9 (pAAV/D377Y-mPCSK9) was delivered via a single tail vein injection⁵². One week later, mice were placed on either a high-cholesterol western diet (Envigo, TD.88137; 42% calories from fat, 0.2% cholesterol) or a matched control diet for 13 weeks. At endpoint, adoptive bone marrow transfer mice were euthanized and Tet2^–/– (CD45.2⁺) or wild-type (CD45.1⁺) neutrophils were isolated from peripheral blood and bone marrow by cell sorting. Cells were processed for scRNA-seq as described below.

Sample preparation and flow cytometry-assisted cell sorting

Mice were euthanized and Blood was taken through cardiac puncture with a 1 ml syringe attached to a 26G needle filled with 50 ml of 0.5 M EDTA. After blood collection, mice were perfused via the right ventricle with 10 ml of PBS to remove circulating blood cells. Tissues, including lung, tumours, muscle, heart, placenta and pancreas, were collected, cut into small pieces, and digested with Liberase TM (Sigma) and DNase I (Sigma) for 30 min at 37 °C. Following digestion, tissues were passed through 70-µm nylon mesh sieves using syringe plungers to obtain single-cell suspensions.

Bone marrow cells were obtained by flushing femurs with PBS containing 2 mM EDTA and 2% FBS using a 23G needle. Spleens were mechanically dissociated through 70-µm mesh filters. For colon isolation, intestines were cleaned, cut longitudinally and washed in PBS. After a 30-min incubation in 100 mM EDTA at 37 °C with shaking to remove epithelial cells, colon tissue was cut and digested in Liberase TM and DNase I for 30 min at 37 °C. Ear (skin) samples were processed by separating the dorsal and ventral sides, cutting them into small pieces, and digesting them in Liberase TM and DNase I for 90 min at 37 °C. The resulting suspensions were filtered as above. Peritoneal lavage was performed by injecting 10 ml of cold PBS into the peritoneal cavity, followed by gentle massage of the abdomen and careful aspiration of the fluid using a syringe and needle. The collected fluid was centrifuged, and the pellet was resuspended in fluorescence-activated cell sorting (FACS) buffer for staining.

For meninges isolation, mice were euthanized and decapitated. The skull was opened along the sagittal midline, and the brain was removed to expose the dura. The meninges were peeled off using fine forceps and placed directly into digestion solution on ice. For the brain, infarcted regions were dissected and digested for 30 min at 37 °C in an enzyme cocktail containing: 50 U ml⁻¹ collagenase; 8.5 U ml⁻¹ dispase; 100 µg ml⁻¹ Nα-tosyl-L-lysine chloromethyl ketone hydrochloride; 5 U ml⁻¹ DNase I in 9.64 ml HBSS without calcium, magnesium, or phenol red (Fisher Scientific, 14175-095). After digestion, brains were ground using a 2 ml glass-glass grinder, filtered through a 70-µm filter, and centrifuged. Cell pellets were resuspended in 7 ml of 35% Percoll, overlaid with 5 ml HBSS to form a gradient, and centrifuged at 800g for 30 min at 4 °C (no brake). The myelin layer and supernatant were discarded, and the final cell pellet was washed and resuspended in FACS buffer.

All single-cell suspensions were lysed in RBC lysis buffer (eBioscience) for 4 min and stained with the following antibodies: CD11b-PE (Clone M1/70, BioLegend, 1:200); LY6G-AF647 (Clone 1A8, eBioscience, 1:200); DAPI (1:10,000). Neutrophils were sorted as live (DAPI-negative), CD11b⁺LY6G⁺ cells using a FACS Aria sorter (BD Biosciences) at the Centro Nacional de Investigaciones Cardiovasculares (CNIC) Cytometry Unit. Bone marrow neutrophils were captured as lineage-negative (B220, CD18, NK.1.1, Ter119, CD3).

Cancer cell culture

The C57BL6/J syngeneic mouse LLC, E0771 breast luminal B and the BALB/c-derived 4T1 breast tumour cell cancer cell lines were from the American Type Culture Collection. The pancreatic adenocarcinoma FC1242 cell line (gift from D. Engle) was derived from Pdx1^cre; Kras^G12D/+; null/+ (KPC) mice. B16-OVAGFP cells were provided by the laboratory of D. Sancho. All cells were cultured in DMEM (Thermofisher) supplemented with 10% FBS (Thermofisher) and 100 μg ml⁻¹ penicillin/streptomycin (Thomas Sci).

In vitro mouse neutrophil culture and analysis

Primary mouse neutrophils were obtained from the femurs and tibias of healthy C57BL/6J mice, or indicated genetically modified mouse model, by centrifugation. Erythrocytes were lysed using Red Blood Cell Lysis Solution (Qiagen; 79217). Cell strainer-filtered single-cell solutions were sorted in BD Aria Cell Sorter as DAPI⁻ CD45⁺CD11b⁺LY6G⁺CD101⁺ mature, and DAPI⁻ CD45⁺CD11b⁺LY6G⁺CD101⁻ immature neutrophils. Cells were seeded in 96-well plates, 50,000 cells per well, and cultured with complemented DMEM medium (vehicle), or with the indicated treatments. G-CSF (574606, BioLegend), TGFβ (7666-MB-005/CF, R&D), IFNβ (581302, BioLegend), CXCL12 (578704, BioLegend), IL-1β (401-ML, R&D Systems), and GM-CSF (315-03-20UG, Thermofisher) were used at a concentration of 10 ng ml⁻¹. Conditioned medium of LLC or KP-PDAC cells was obtained after 24 h culture of 80% confluence cells. Neutrophils were collected after 24 h or 48 h of treatment, and flow cytometry was performed using the following antibodies, all diluted 1:200 unless indicated otherwise: CCR5-BUV615-P (BD Biosciences, 752321, clone C34-3448), CD101-PE-Cy7 (eBioscience, 25-1011-82, clone MOUSHI 101), CD106-BUV563 (BD Biosciences, 741246, clone 429), CD115-BUV737 (BD Biosciences, 750948, clone AFS98), CD11b-BV510 (BioLegend, 101263, clone M1/70), CD11b-PE (BioLegend, 101208, clone M1/70), CD14-APC-Cy7 (BioLegend, 123318, clone Sa14-2), CD150-PE-Cy5 (BioLegend, 115911, clone TC15-12F12.2), CD16/32-PerCP-Cy5.5 (BioLegend, 101324, clone 93), CD274-BV421 (BioLegend, 124315, clone 10F-9G2), CD44-BV570 (BioLegend, 103037, clone IM7), CD45-APC (BioLegend, 103112, clone 30F11), CD74-BUV661 (BD Biosciences, 741572, clone In-1), KIT-BV605 (BioLegend, 135121, clone ACK2), CX3CR1-FITC (BioLegend, 149020, clone SA011F11), DC-Trail-R1-biotinylated (R&D Systems, BAF2378, polyclonal), I-A/I-E-BUV496 (BD Biosciences, 750281, clone M5/114.15.2), ICAM1-PE-Dazzle 594 (BioLegend, 1161130, clone YN1/1.7.4), LY6C-BV711 (BioLegend, 128037, clone HK1.4), LY6G-PE (BioLegend, 127608, clone 1A8), Sca1-BUV395 (BD Biosciences, 563990, clone D7), TLR4-BV786 (BD Biosciences, 741015, clone MTS510). Streptavidin-BV650 (BioLegend, 405231) was included at 1:500.

Secondary staining was performed with Streptavidin-BV650 (Biolegend). Cells were analysed in a SymphonyA4 Flow Cytometer. The data were analysed using FlowJo v.10 software. FlowAI⁵³ was used for quality control of flow data, followed by dimensionality reduction using the UMAP_R plugin⁵⁴. Initial clusterization was performed with FlowSOM⁵⁵ and ClusterExplorerPlugin, and UMAP parameters were embedded in each sample for statistical analysis of neutrophil phenotypes.

RNA isolation, reverse transcription PCR

Total RNA was prepared with the RNA Extraction RNeasy Plus Mini-kit (QIAGEN) and RNA was reverse-transcribed with the High-Capacity cDNA Reverse Transcription kit (Applied Biosystems) according to the manufacturer’s protocol. Real-time quantitative PCR (SYBR-green, Applied Biosystems) assays were performed with an Applied Biosystems 7900HT Fast Real-Time PCR System sequencer detector. Expression was normalized to the expression of the 36b4 housekeeping gene. Primer sequences are as follows: 36b4: forward 5′-ACTGGTCTAGGACCCGAGAAG-3′, reverse 5′- TCCCACCTTGTCTCCAGTCT-3′; Ptgs2: forward 5′-TGAGCAACTATTCCAAACCAGC-3′, reverse 5′-GCACGTAGTCTTCGATCACTATC-3′; Nr4a1: forward 5′-TTGAGTTCGGCAAGCCTACC-3′, reverse 5′-GTGTACCCGTCCATGAAGGTG-3′; Il1b: forward 5′-AGTGAGGAGAATGACCTGTTC-3′, reverse 5′-CGAGATGCTGCTGTGAGATT-3′; Tnfaip3: forward 5′-GAACAGCGATCAGGCCAGG-3′, reverse 5′-GGACAGTTGGGTGTCTCACATT-3′; Cebpe: forward 5′-GCAGCCACTTGAGTTCTCAGG-3′, reverse 5′GATGTAGGCGGAGAGGTCGAT-3′; Ltf: forward 5′-TGAGGCCCTTGGACTCTGT-3′, reverse 5′-ACCCACTTTTCTCATCTCGTTC-3′.

Functional assays

T cell cytotoxicity assay

B16F10–OVA-GFP (10⁴ cells) were seeded in 96-well culture dishes for 24 h, in RMPI medium with glutamine (Thermofisher) containing 10% heat-inactivated FBS (Thermofisher), 100 μg ml⁻¹ penicillin/streptomycin (Thomas Sci); 200 nM glutamine, 1% non-essential amino acids (MEM amino acids; Gibco), 1% sodium pyruvate (Gibco) and 0.01% β-mercaptoethanol (Gibco). Neutrophils from sorting or in vitro cultures were co-culture at a 2:1 ratio with SIINFEKL-activated OT-1 T effector cells for 3 h. Neutrophil-OT-I cells were then seeded on top of B16F10-OVAGFP target cells 1:2 ratio. After 24 h, live cells were stained with 0.4 g l⁻¹ crystal violet (Sigma-Aldrich, HT90132). The area covered by target cells was quantified from micrographs of the plates using the ImageJ software.

In vivo Matrigel plug assay

Fifty thousand neutrophils sorted from tissues of interest were resuspended in 500 µl of growth-factor-depleted Matrigel (Corning) and injected subcutaneously in the lower back of anaesthetized mice to form plugs. At days 3 and 7 after implantation, the same number of sorted neutrophils was resuspended in 50 µl of PBS and injected directly into the plug respectively. On day 9 after implantation, Doppler laser perfusion imaging was performed at the lower back region that contained the Matrigel plugs. One ROI was defined for each observable Matrigel plug, and the amount of flux variation in each ROI was quantified. Only ROIs that were not obscured by hair regrowth were used.

Chemotaxis assay

Chemotaxis assays were performed as described⁴⁷. In brief, bone marrow neutrophils were plated in 6.5 mm polycarbonate transwells with 5-mm pores (Corning) in RPMI medium 48 h after cytokine treatment. 20 ng ml⁻¹ CXCL1 (R&D) was added to the bottom well. Transwells were incubated for 2 h at 37 °C and transmigrated cells were collected from the bottom well and stained for cytometric analysis. The number of transmigrated cells was assessed by the presence of a known number of Truecount beads (BD Biosciences).

3D Doppler imaging of tumour vascularization

Subcutaneous LLC tumour vascularization was imaged using Vevo Imaging Systems once they reached 500 mm³. In brief, mice were anaesthetized in an isoflurane vapourizer chamber, and the backs were thoroughly shaved. The mice were placed in the imaging platform and images were captured using the power colour Doppler-3D mode. A total of 100 images were captured to generate a 3D reconstruction of the vasculature. Vevo LAB software was used to calculate the Volume and per cent vascularization of tumours. Per cent vascularization is determined by calculating the percentage of pixels in the volume that have a power Doppler signal associated with them, the presence of this signal indicates the presence of blood flow.

NET formation assay

Forty-eight hours after cytokine treatment, 5 × 10⁴ bone marrow neutrophils were plated with RPMI medium on poly-l-lysine-covered 8-well μ-Slides (Ibidi), and left 30 min to adhere. Subsequently, cells were incubated for 2 h with 100 nM PMA or vehicle. Cells were then fixed using 4% PFA for 10 min, permeabilized with PBS with 0.1% Triton X-100, 1% goat serum plus 5% BSA and stained with antibodies against cit-H3, DNA (Sytox-green, Molecular Probes) and MPO. Whole-slide z-stack tilescan images were acquired with a Leica SP5 confocal microscope, and analysed using Imaris software (v.9.5, Bitplane)⁴⁷.

Bacterial killing assay

Forty-eight hours after cytokine treatment, bone marrow neutrophils were resuspended in fresh medium along with live S. aureus (ATCC) that were grown in tryptic soy broth. For the in vitro assays, neutrophils and bacteria (10⁴ CFU in 200 μl) were incubated at 37 °C for 60 min. The cells are then plated onto tryptic soy plates in a serial dilution. Bacterial colonies on the plates were counted the following day.

Phagocytosis assay

Forty-eight hours after cytokine treatment, bone marrow neutrophils were resuspended in fresh medium along with fluorescent latex beads (SIGMA) followed by flow cytometric analyses.

Analysis of human neutrophils

Isolation and expansion of human bone marrow CD34⁺ HSPCs

Bone marrow samples were obtained from healthy donors under informed consent approved by the ethics committee of the University Hospital Tübingen. CD34⁺ haematopoietic stem and progenitor cells (HSPCs) were isolated through Ficoll gradient centrifugation followed by magnetic bead-based separation using the EasySep Human CD34⁺ Cell Selection Kit II (Stem Cell Technologies, 17856). CD34⁺ cells (n = 4; purity 95.4 ± 1.9%) were cultured at a density of 5 × 10⁵ cells per ml in StemSpan SFEM II haematopoietic stem cell medium (Stem Cell Technologies, 09655), supplemented with 1% penicillin/streptomycin, 20 ng ml⁻¹ IL-3, 20 ng ml⁻¹ IL-6, 20 ng ml⁻¹ TPO, 50 ng ml⁻¹ SCF and 50 ng ml⁻¹ FLT-3L (all cytokines purchased from R&D Systems). Cells were cultured under standard conditions (37 °C, 5% CO₂) and frozen for future use.

For granulocytic differentiation in vitro, cells were seeded at a density of 2 × 10⁵ cells per ml. During the first 8 days of differentiation (days 0–7), cells were maintained in a myeloid cell expansion medium—RPMI 1640 supplemented with 10% FCS, 1% penicillin/streptomycin, 5 ng ml⁻¹ SCF, 5 ng ml⁻¹ IL-3 and 1 ng ml⁻¹ G-CSF. The medium was changed every two days. On day 8 of culture, the medium was replaced with a granulocytic cell differentiation medium (RPMI 1640 supplemented with 10% FCS, 1% penicillin/streptomycin and 1 ng ml⁻¹ G-CSF). The medium was changed every 2 days until day 14. On day 13 of differentiation, cells were collected and counted. 800,000 cells were lysed for RNA isolation, 50,000 cells for FACS, and 40,000 cells for cytospins. The remaining cells were resuspended in fresh granulocytic differentiation medium at a seed density of 2 × 10⁵ cells per ml and divided into 4 groups. Group one was maintained in granulocytic differentiation medium, group two was treated with 10 ng ml⁻¹ TGFβ, group three was treated with 10 ng ml⁻¹ IFNβ (refreshed after 24 h), and group 4 was treated with 10 ng ml⁻¹ GM-CSF. RNA-seq analyses were performed 48 h after stimulation.

HOXB8 cell cultures and differentiation

HOXB8-immortalized myeloid progenitors were routinely tested for mycoplasma contamination and cultured in RPMI 1640 medium supplemented with 10% fetal calf serum (FCS), 10 μM β-mercaptoethanol (Thermo Fisher Scientific), 4% supernatant from SCF-producing CHO cells, 1% penicillin/streptomycin and 1 μM β-oestradiol (Sigma-Aldrich) to maintain the progenitor state. Neutrophil differentiation was initiated by β-oestradiol withdrawal and continued culture in medium supplemented with 1% SCF-containing supernatant. Differentiation into neutrophils was achieved by culturing cells in RPMI 1640 medium containing 10% FCS, 30 μM β-mercaptoethanol, 4% SCF supernatant, and 20 ng ml⁻¹ granulocyte colony-stimulating factor (G-CSF) under standard tissue culture conditions (37 °C, 5% CO₂).

CRISPR–Cas9-mediated knockout

Knockouts of selected transcription factors in HOXB8 progenitors have been previously described²¹. In brief, HOXB8 progenitors were transduced with lentiCas9-v2 lentiviral vectors encoding guide RNAs (gRNAs) targeting the following genes and exons: Cebpb (exon 1; gRNA: AGGCTCACGTAACCGTAGT); Klf6 (exon 1; gRNA: TCGCTGTCGGGAAAACAGGG); Runx1 (exon 3; gRNA: TAGCGAGATTCAACGACCTC); Rfx2 (exon 5; gRNA: CTGCTGGGGGCGTAAAGCTG); Relb (exon 4; gRNA: CTGCACGGACGGCGTCTGCA); Irf5 (exon 2; gRNA: ACCCTGGCGCCATGCCACGAGG); and Junb (exon 1; gRNA: GGAACCGCAGACCGTACCGG).

JUNB overexpression

Lentiviral vectors for JUNB overexpression were generated by transient transfection of HEK293T cells using the calcium phosphate precipitation method. Cells were co-transfected with: (1) a transfer plasmid containing the Junb cDNA under the control of the human PGK promoter; (2) packaging plasmid psPax2; and (3) envelope plasmid pMD2.G encoding VSV-G. The medium was replaced 24 h after transfection. At 72 h, virus-containing supernatant was collected, clarified by centrifugation (2,000 rpm, 5 min, 4 °C), filtered (0.45 μm), and concentrated via ultracentrifugation (26,000 rpm, 2 h, 4 °C). Viral pellets were resuspended in cold PBS, aliquoted, and stored at −80 °C.

Lentiviral transduction of HOXB8

HOXB8 progenitors were transduced by spinoculation. In brief, 5 × 10⁵ cells were plated per well in 6-well plates with 1 ml of medium. Lentiviral particles were added at a multiplicity of infection (MOI) of 11.24 for the vector pRRL-hPGK-JUNB-IRES-eGFP and MOI = 1.8 for the pRRL-hPGK-IRES-eGFP empty vector, and cells were centrifuged at 1,000g for 90 min at 30 °C. Following transduction, cells were collected, washed, and resuspended in fresh culture medium at a final concentration of 5 × 10⁴ cells per ml.

Plasmid construction

To construct the JUNB expression vector, the human PGK (hPGK) promoter was PCR-amplified with ClaI and XbaI restriction sites and cloned into the ClaI/XbaI sites of the pRRL-CMV-IRES-eGFP vector, replacing the CMV promoter. The Junb coding sequence was amplified from mouse cDNA using primers containing BglII and XhoI sites and inserted into the BamHI and XhoI sites of the modified vector. Cloning was performed using the following primers: hPGK forward (ClaI): 5′-TTTTTTATCGATGGGTAGGGGAGGCGCTTT-3′; hPGK reverse (XbaI): 5′-TTTTTTTTAGACGAAAGGCCCGGAGATGA-3′; Junb forward (BglII): 5′-TTTTTTAGATCTGCCACCATGTGCACGAAAATGGAACA-3′; Junbreverse (XhoI): 5′-TTTTTTCTCGAGTCAGAAGGCGTGTCCCTT-3.

Culture of HOXB8 cells

For the flow cytometry and bulk RNA-sequencing experiments, HOXB8 progenitors at day 3.5 of differentiation were seeded in 96-well plates at a density of 50,000 cells per well. Cells were cultured in complete RPMI medium (vehicle) or treated with GM-CSF (10 ng ml⁻¹) for 48 h, following the same conditions described for primary bone marrow neutrophil cultures. Vehicle or GM-CSF treated cells were collected at 48 h after treatment for the analysis.

Bulk RNA sequencing of mice and human-derived neutrophils

RNA from isolated mouse neutrophils was extracted using RNAeasy micro kit (Quiagen). RNA quality was checked using capillary electrophoresis (Agilent). Samples were submitted for whole RNA next generation sequencing in the Genomics Unit of CNIC. Total RNA (200 ng) was used to generate barcoded RNA-sequencing libraries using the NEBNext Ultra RNA Library preparation kit (New England Biolabs). Libraries were sequenced with HiSeq2500 (Illumina) to generate 50-nucleotide single reads, with a minimum of 8 million reads per sample. For RNA-seq of human-derived neutrophils, we isolated RNA from a total of 800,000 differentiated neutrophils collected on day 13 and 15. We used the NucleoSpin RNA Mini Kit (Macherey-Nagel, 740955.50), following the manufacturer’s instructions. RNA concentration was assessed with Qubit 2.0 (Thermo Fisher), and a total of 400 ng RNA was sequenced. RNA integrity was assessed using Agilent Bioanalyzer 2100, with RNA integrity number (RIN) between 9.8 and 10.0. RNA samples were processed by Novogene for library preparation and sequencing, and all samples passed the quality control criteria. Strand-specific libraries were generated on the basis of Novogene’s standard protocol. Samples were sequenced on an Illumina platform to produce 150 bp pairwise reads (PE150) per sample.

FastQ files for each sample were obtained using CASAVA (v.1.8) software (Illumina). Reads were further processed with RTA (v.1.18.66.3). FastQ files for each sample were obtained using bcl2fastq (v.2.20.0.422) software (Illumina). Sequencing reads were further processed as follows: Illumina adapters were trimmed and low-quality reads removed with Cutadapt (v.4.9)⁵⁶ (mismatch rate = 1 mismatch every 10 bp, overlap = 5 bp, minimum read length = 30 bp). Quality control of the processed reads was done with fastQC (v.0.12.1). RSEM (v.1.3.1) was used to quantify expression levels against the mouse genome reference GRCm38 or the human genome reference GRCh38, depending on the analysis⁵⁷ (default options). The processing of the counts and differential expression analysis was performed using limma (v.3.32.2)⁵⁸ and EdgeR (v.3.20.1)⁵⁹) which were also used to perform pairwise differential expression analyses. To identify genes whose expression significantly varies across conditions, we applied a Likelihood Ratio Test (LRT) using DESeq2 (v.1.30.1)⁶⁰, allowing the detection of global effects of a factor without the need to specify individual contrasts. The resulting significant genes were then clustered using the k-means algorithm from the stats package (v.4.0.3).

Single-cell transcriptomics on mouse neutrophils

scRNA-seq of sorted tissue neutrophils

For single-cell analysis, all samples were collected between ZT1 and ZT5. Tissues were dissected and dissociated into a single-cell suspension by enzyme digestion. The resulting suspensions were filtered through cell strainers, and sorted in BD Aria Cell Sorter as DAPI⁻CD11b⁺LY6G⁺ cells, and loaded into a BD Rhapsody cartridge. For the generation of single-cell whole-transcriptomes, we used a BD Rhapsody system according to the manufacturer’s instructions. In brief, cell suspensions from each condition were incubated with Sample Tags (BD) for 20 min at room temperature. Cells were then washed three times and pooled in a single tube. Cell viability and concentration were assessed using a Countess III cell counter (Thermo Fisher). Sixty thousand cells were loaded into a Rhapsody Single Cell Analysis System cartridge. Cell capture and cDNA synthesis were performed according to manufacturer’s instructions; cells were isolated into nanowells by gravity, then cells were lysed and mRNAs together with sample tags oligonucleotides were released and captured by the beads present in the nanowells. Each bead contained a unique oligonucleotide named ‘cell label’ to identify each individual bead. All beads present in the cartridge were collected and cDNA synthesis took place in a single reaction. At this point, each cDNA and Sample Tag oligonucleotide were attached to its corresponding cell label oligonucleotide. Two separated indexed libraries were prepared for whole-transcriptome analysis and sample tag demultiplexing following the manufacturer’s instructions. The average size of the libraries was calculated using the 2100 Bioanalyzer (Agilent), and the concentration was determined using the Qubit fluorometer (Thermofisher). Finally, libraries were combined and sequenced together in a paired-end run (60 × 42) using a NextSeq 2000 system (Illumina) and a P2 flow cell. Output files were processed with NextSeq 1000/2000 Control Software Suite v.1.4.1. FastQ files for each sample were obtained using BCL Convert v.3.6.3 software (Illumina).

Construction of NeuMap and projection of external data

Rhapsody analysis pipeline v.1.9.1 was run locally. This pipeline includes steps for, alignment to mouse genome reference (GRCm38 with the gencodevM19-20181206) quantification and filtering of low-quality cells and tagging of doublets, which were also filtered out of the downstream analyses⁶¹. After BD Rhapsody’s pipeline automatic quality filtering, a second round was performed manually, where cells with a mitochondrial content over 20% or with over 300 total gene counts were discarded. Cell Annotation was performed using R package SingleR and the Immgen database for each dataset individually. All subsequent downstream analyses were implemented using R (v.4.0.3) and the package Seurat (v.4.0.5)⁶¹. The Seurat suite was used to integrate the neutrophils from all datasets using Seurat’s integration implementation. This method uses common sources of variation across the different datasets and aligns the cells so those in similar biological states cluster together. The integrated dataset was used to perform the unbiased cluster analysis and the construction of NeuMap. Additionally, we used the integrated NeuMap to generate a reference which we later used to analyse additional and external datasets by projecting cells onto our reference and annotate the new data using our custom labels using Seurat’s MapQuery and TransferData. This method is technology agnostic, so we could reliably project cells from external datasets sequenced in different platforms onto NeuMap^62,63.

Definition of hubs

Functional hubs were selected by performing unbiased clustering at different resolutions using Seurat’s function FindClusters(). Resolutions used ranged between 0.05 and 0.3. Clusters from different resolutions were selected because they best represented the expression of functional signatures projected onto NeuMap. Areas shown in the figures correspond to the q15 quantile of the KMASS algorithm, which calculates the density of cells in specific areas. For clarity, hubs in figures are shown as the area with the accumulation of 85% of cells for each selected cluster/hub. Analyses were performed on the complete set of cells for each cluster or hub. The FindAllMarkers() function from Seurat was used to calculate DEGs across the hubs. Only genes detected in a minimum of 25% of the cells and with an average of at least 0.25-fold difference (log scale) between the groups in either of the groups were tested.

Kernel density estimation

The MASS R package (v.7.3.61) was used for two-dimensional kernel density estimation (K-mass score), with n = 100 grid points in each direction.

Signature projection

The signatures used for illustration of functional states are contained in Supplementary Table 2. All signatures were calculated by Seurat’s AddModuleScore() function. We used two different sources for the functional signatures: (1) previous publications, for which we provide the whole list of genes reported and used in Supplementary Table 2; and (2) public databases such as gene ontology (GO) and gene set enrichment analysis (GSEA). In those cases the whole gene list from the functional category was tested. For signatures from human data, human genes were mapped to their corresponding mouse homologue to calculate the enrichment score using the Mouse Genome Informatics (MGI) database. We used Seurat R package (v.5.2.1) AddModuleScore() function to calculate the scores. To generate visualization heat maps across NeuMap hubs, we first calculated enrichment scores for each cell. Scores were then averaged by hub and scaled per signature for comparison. To assess whether gene signature scores significantly differed across NeuMap hubs, we applied the Kruskal–Wallis test to each signature, testing the null hypothesis that score distributions were identical across hubs. The resulting test statistics were compared to a chi-squared distribution, with degrees of freedom equal to the number of hubs minus one. To correct for multiple comparisons, we adjusted P values using the Benjamini–Hochberg false discovery rate method.

Velocyto analysis

The analysis of expression dynamics in scRNA-seq data was performed using velocyto (v.0.17.17)²⁴, a package that allows estimating RNA velocities distinguishing between spliced and unspliced mRNAs in standard scRNA-seq protocols. The command line tool in Python implementation was adapted to be able to work with BAM files generated by BD Rhapsody, using samtools⁶⁴ to format the files, mainly by removing all possible alignments with antibodies and renaming the UMI barcode tag to ‘UB’ instead of ‘MA’. Velocyto was then executed with default parameters and the GRCm38 reference genome with the gencodevM19-20181206 transcriptome annotation. After concatenation of the spliced and unspliced data from all experiments, the results were merged with the outputs from single-cell analyses performed with Seurat in R, and scVelo⁶⁵ was used for further processing. Pre-processing included gene selection by detection (the minimum number of both unspliced and spliced counts was set to 30), and by variability (keep 2,000 highly variable genes (HVGs)), normalization, and log1p transformation. First and second order moments were computed among 50 nearest neighbours in the principal component analysis (PCA) space using 30 components. Cell-based RNA velocities were estimated by modelling the transcriptional dynamics of splicing kinetics using the stochastic model available in scVelo. Finally, these velocities were projected onto the previously computed UMAP and visualized at the cellular level or as velocity vector fields through streamlines.

In some experiments we performed Pseudotime analysis. Samples were pre-processed using the standard Monocle3 pipeline. To address batch effects, samples were integrated using the Batchelor algorithm (v.1.20.0). Dimensionality reduction and clustering were performed within Monocle3 (v.1.3.7), and pseudotime values were computed for the integrated dataset. To evaluate the significance of differences in pseudotime values between Cre⁻ control and Tgfbr2-mutant immature cells, we applied a non-parametric Wilcoxon rank-sum test with continuity correction (W = 373,675; P value = 0.005548).

Comparison of cell distribution in the different hubs

To test differences in how neutrophils are distributed in the different hubs, we classified cell hubs for both mutant and control samples on the basis of the k-nearest neighbours algorithm (k = 5) of the cells projected onto Neumap. The observed hub proportions were calculated for each sample, and differences were determined by subtracting the proportions in control from the proportions in mutant cells. To assess the statistical significance of these differences, we used a bootstrap approach. For this, we generated a null distribution of hub proportion differences by merging each control–mutant pair into a mixed population. From this combined dataset, 10,000 resampled pairs were drawn with replacement, matching the sample sizes of the original control and mutant datasets. The differences in hub proportions between the resampled mutant and control groups were then calculated. The null distributions for each hub were verified to have a mean of 0.0, as expected under the null hypothesis. Finally, P values for the observed differences were computed by determining the fraction of resampled differences that were as extreme as or more extreme than the observed differences. To estimate 95% confidence intervals, the quantiles corresponding to the 0.025 and 0.975 percentiles of the null distributions were calculated.

Single-cell multiome using Dogma-seq

To simultaneously profile chromatin accessibility and gene expression at single-cell resolution, we used the Chromium Next GEM Single Cell Multiome ATAC + Gene Expression platform (10x Genomics). We collected neutrophils from bone marrow, spleen and lung (dataset 1) and LLC and spleen with LPS (dataset 2) from 8-to-12-week-old C57BL/6 mouse healthy blood was sequenced in both datasets as a quality control reference. Single-cell suspensions were prepared as described above. After staining, cells were washed, resuspended in sorting buffer, and incubated with DAPI (NBP2311561, Novus Biologicals) for 15 min prior to sorting.

Live CD11b⁺LY6G⁺ neutrophils were sorted in equal proportions from each organ. Cells were then pooled and lysed in 100 µl of cold DIG lysis buffer (20 mM Tris-HCl, 150 mM NaCl, 3 mM MgCl₂, 0.005% digitonin, 2 U μl⁻¹ RNase inhibitor) for 5 min on ice. Lysis was quenched with 1 ml of cold DIG wash buffer, followed by centrifugation at 500g for 5 min. Nuclei were resuspended in 100 µl of 10x Genomics Nuclei Buffer supplemented with 1 mM DTT (Sigma) and 2 U µl⁻¹ RNase inhibitor (Roche) to a final concentration of 3,400 nuclei per µl. After additional washes and centrifugation, samples were processed for library preparation at the Yale Center for Genome Analysis.

Library construction was performed following the 10x Genomics protocol (Chromium Next GEM Multiome ATAC + GEX v.1.1, CG000338 rev. E). In brief, nuclei underwent transposition using the ATAC transposition mix and were loaded onto the Chromium Controller for GEM generation, barcoding, and reverse transcription. Separate libraries were constructed for ATAC and gene expression using standard amplification and indexing steps. Libraries were quantified using Bioanalyzer (Agilent) and Qubit (ThermoFisher), pooled, and sequenced on an Illumina NovaSeq 6000 platform (paired-end, 150 bp) with a target depth of 75 million reads per sample.

Data processing

Initial quality control and cell filtering

DNA accessibility and gene expression from each cell were analysed simultaneously using Seurat (v.4.0.5 and v.4.3)⁶¹ and Signac (1.14.0)⁶⁶ R packages. Per cell quality control metrics were evaluated using the DNA accessibility, and transcriptional data were obtained. Cells that did not pass the following criteria were removed from downstream analysis: number of counts in the ATAC data 100 < (nCount_ATAC) < 100,000; number of counts in the gene expression data 100 < (nCount_RNA) < 5,000; number of genes in the gene expression data 100 < (nFeature_RNA) < 2,500; ratio of mono-nucleosomal to nucleosome-free fragments (nucleosome_signal) < 2; ratio of fragments centred at the transcription start site (TSS) to fragments in TSS-flanking regions > 1; percentage of mitochondrial gene expression < 5.

ATAC data annotation

R packages Signac (v.1.14.0) and Seurat (v.5.1.0) were used to analyse single-cell chromatin data and gene expression, respectively. Full genome sequences for M. musculus (mouse) were used as provided by UCSC (mm10, based on GRCm38.p6), and annotated using Ensembl M. musculusannotations v.79.

Cell-type identification and neutrophil subset classification

We used R package SingleR (v.2.8.0)⁶⁷ to annotate cell types against the ImmGen database⁶⁸. Cells annotated as ‘neutrophils’ were subset and re-analysed by running a new round of FindVariableFeatures() in which outlier features were identified and ScaleData() to re-scale the expression of the neutrophil subset.

Mapping onto NeuMap

We used Seurat v.5.1. FindTransferAnchors() function to identify pairwise correspondences (anchors) between the reference and query datasets using the transcriptomics data. This function uses canonical correlation analysis and mutual nearest neighbours to identify cells with similar gene expression profiles across the two datasets. The top 2,000 variable features shared between the reference and query datasets were used for anchor identification.

The query dataset was mapped onto the NeuMap reference using the MapQuery() function. This step projected the query cells from the Dogma-seq into NeuMap embedding space, allowing direct comparison and visualization of the dogma cells relative to NeuMap. Hub annotations from NeuMap were transferred to the query dataset using the TransferData() function. This function predicts cell labels for each query cell on the basis of the similarity scores computed from the anchors. The predicted labels were assigned to the query dataset, enabling downstream analysis of chromatin state in cells from each hub. Additionally, we assessed the confidence scores provided by TransferData() for each predicted label, retaining only high-confidence predictions (predicted.id.score ≥ 0.7) for downstream analysis.

Merging of the datasets and peak calling

We created a common peak set, and quantified this peak set in each experiment using Signac (v.1.14.0) and GenomicRanges (v.1.58.0)⁶⁹ prior to merging the objects. Once both datasets contained an assay with the same set of features, we used Seurat (v.5.2.1) R package to merge the datasets.

We used the Signac R package (v.1.14.0) to call peaks with the CallPeaks() function. The CallPeaks() function used MACS2 (v.2.2.9.1)⁷⁰ to run. Peaks were called for cells assigned to each hub separately. Only cells with a predicted.id.score ≥ 0.7 were retained for peak calling. Peaks on nonstandard chromosomes and in genomic blacklist regions were removed. After quality control and predicted score filtering 1,962 for dataset 1 and 8,155 neutrophils remained for dataset 2.

ATAC data processing

R package Signac (v.1.14.0) standard processing pipeline was applied to the combined data: term frequency-inverse document frequency normalization was applied via RunTFIDF(), top features were identified using FindTopFeatures() with a minimum cut-off of 5 and singular value decomposition was performed on the normalized data running RunSVD().

Differential peak analysis

Differential accessibility peaks were identified using FindAllMarkers(), considering only positive markers and a minimum percentage of cells expressing the feature (min.pct = 0.1). The closest genes to the differentially accessible peaks were annotated using the ClosestFeature(). The results were merged and filtered to retain significant peaks and marked for uniqueness.

Motif enrichment analysis

A position frequency matrix set was retrieved from the JASPAR2020 database via the homonim R package (v.0.99.10), filtering for vertebrate transcription factors in the CORE collection. Motif information was added to the dataset with the BSgenome.Mmusculus.UCSC.mm10 genome as reference. Enriched motifs in the differentially accessible peaks per hub were then identified.

Transcription factor activity

Chromatin accessibility variability analysis was performed using R package chromVAR (v.1.28.0)⁷¹, with the BSgenome.Mmusculus.UCSC.mm10 genome as the reference. This step computes motif activity scores for each cell, representing the inferred transcription factor activity based on chromatin accessibility.

Single-cell transcriptomics on human neutrophils

Samples collection and processing

Human samples were collected in Renji Hospital, Shanghai, China, under the Renji Hospital Ethics Committee protocol KY2024-090-B, in accordance with the Declaration of Helsinki, following informed consent from all participants. Samples were collected from healthy donors, patients, or perfused organ donors. Specifically, healthy donor samples (bone marrow, peripheral blood and umbilical cord blood) were randomly collected without self-selection or recruitment bias. Other healthy tissues were obtained from anonymous acute-death donors without chronic inflammation to minimize the confounding effects of death shock on the organs. Systemic lupus erythematosus (SLE) patient samples (umbilical cord and peripheral blood) were randomly selected from pregnant patients with an active disease state (SLEPDAI > 5) and without other chronic inflammatory comorbidities.

Blood and bone marrow were collected in BD vacutainer K2E (EDTA) tubes (BD Healthcare, 367525) to prevent coagulation. Erythrocytes were lysed in 5–10 ml 1× red blood cells (RBC) lysis buffer (diluted from 10× BD Pharm Lyse, 555899) for 5 min for twice to deplete erythrocytes and then washed and re-suspension. Spleen, lung, omentum, mesentary fat, perirenal fat, liver, colon and rectum tissues were minced into small pieces and digested for 30 min at 37 °C in a mixture of collagenase IV (385 U ml⁻¹, Sigma) and DNase I (2.5 mg ml⁻¹, Sigma) and the samples were homogenized into single-cell suspension using syringe plungers and passed through 70-μm cell strainers (15-1070, BIOLOGIX). Then the samples were lysed in 2 ml 1× RBC lysis buffer (diluted from 10× BD Pharm Lyse, 555899) for 3 min to deplete erythrocytes and then washed and resuspended. Endometrium was cut into small pieces and enzymatically digested with the Tumor Dissociation Kit (130-095-929, Miltenyi Biotec). After digestion, the cell suspension was filtered through 70-μm cell strainers and subjected to a 3-min erythrocyte lysis with 2 ml 1× RBC lysis buffer, followed by washing and re-suspension.All single-cell suspensions were incubated with Fc-blocker (Human TruStain FcXTM, 422302, Biolegend) for 30 min on ice, then stained for 30 min at 4 °C in the dark with Fixable Viability Stain 700 (564997, BD Biosciences) (1:1,000), and the following antibodies: Anti-CD19-PE-Cy7 (clone HIB19, BioLegend, 302216; 1:200); anti-CD3-PE-Cy7 (clone HIT3a, BioLegend, 300316; 1:200); anti-CD45-APC-Cy7 (clone HI30, BioLegend, 304014; 1:200); anti-CD56-PE-Cy7 (clone 5.1H11, BioLegend, 362510; 1:200). All antibodies were used at 1:200 dilution unless otherwise indicated. After washing with FACS buffer, cells were sorted on a FACS Aria III cell sorter (BD Biosciences). After washing with FACS buffer, cells were sorted on a FACS Aria III cell sorter (BD Biosciences).

Library preparation

Single-cell suspensions were processed on the BD Rhapsody Express System (BD Biosciences). In brief, cells and beads were loaded onto the BD Rhapsody cartridge. Lysis, reverse transcription and exonuclease I digestion were performed using BD Rhapsody Enhanced Cartridge Reagent Kit (BD Bioscience, 664887) and the BD Rhapsody cDNA Kit (BD Bioscience, 633773). The whole-transcriptome libraries were prepared by following the BD Rhapsody single-cell whole-transcriptome amplification workflow with the BD Rhapsody WTA Amplification Kit (BD Bioscience, 633801), including random priming and extension (RPE), RPE amplification PCR and whole-transcriptome amplification index PCR. Libraries were quantified using a High Sensitivity DNA chip (Agilent) on a Bioanalyzer 2100 and the Qubit High Sensitivity DNA assay (Thermo Fisher Scientific) and then sequenced on a NovaSeq X Plus (Illumina). Raw sequencing data (.fqstq files) were processed with the BD Rhapsody analysis pipeline.

Data processing and cell annotation

Seurat v.5.2.1 package standard pipeline was used for the analysis of the single-cell data. The percentage of mitochondrial content was calculated and cell cycle gene expression scores were obtaned using the cell-cycle gene list by Tirosh et al.⁷² and DAM stage genes⁷³ expression scores. Cells with a percentage of mitochondrial content over 20% were removed from downstream analysis. Likewise cells with a number of detected features below 300 were removed. Cells were manually annotated by selecting the clusters in each dataset with highest expression score for known neutrophil markers in different states as described^12,74.

Data integration

scRNA-seq datasets were integrated using the reciprocal principal component analysis (RPCA) method implemented in Seurat v.5.2.1 package. Data normalization and identification of variable features (n = 2,000) were performed independently for each dataset using variance stabilizing transformation. Integration anchors were identified using RPCA reduction with k.anchor=15, and datasets were integrated with k.weight=50. Principal component analysis was performed on the integrated data using the top 15 components for UMAP dimensionality reduction (seed = 42).

Gene module scoring and hub identification

We used Seurat v.5.2.1 AddScoreModules() to assess the activity of specific gene sets within cell clusters. Gene lists of interest shown here were obtained from public data and repositories (Supplementary Table 2), from ref. ¹⁸. For each gene list, we used AddModuleScore() to calculate the aggregated module score against a set of control genes with similar expression, thus ensuring the score was not biased by overall expression levels. To identify human to mouse hub similarities gene lists of interest were obtained from NeuMap hub gene lists. For each gene list, we used AddModuleScore() to calculate the aggregated module score against a set of control genes with similar expression. For all gene module scoring, we reduced the control gene parameter to 80 to ensure sufficient background correction while maintaining computational feasibility within the reduced gene number in the integrated dataset.

Clustering was performed using the Leiden algorithm at several resolution. Clusters with high expression of functional scores were selected as hubs, keeping only cells uniquely assigned to each specific cluster or hub. This approach led to the identification of six human NeuMap regions shown in Supplementary Fig. 4d. For visualization purposes, we used two-dimensional kernel density estimation on UMAP coordinates. For figure visualization, we calculated density values per hub using the kde2d function from the MASS R package (n = 100 grid points) (v.7.3.61), selected cells above the 80th percentile of density (top 20% most accumulated cells) to define core regions representative of each defined hub. The FindAllMarkers() function with default parameters from Seurat was used to calculate DEGs shown is Supplementary Table 4 across the human hubs. Spatial boundaries around these high-density regions were computed using the concaveman algorithm (v.1.1.0) to generate concave hull polygons that capture the geometric extent of each cell state cluster. For the stacked bars plots, all samples are downsampled to 1,000 cells for consistency with previous mouse NeuMap bars and comparability among samples with varying sizes.

Spatial transcriptomics using Visium OCT

Visualization of gene expression in naive lung (n = 1), tumour-bearing lungs (n = 2) and flu-infected lung (n = 1) was performed using the 10x Visium Spatial Gene Expression Kit (10x Genomics; PN1000184) following the manufacturer’s protocol. The OCT blocks were cut using a cryostat (Leica; PN-CM1520) and first cuts were used for RNA extraction (Qiagen; PN-74034), and RNA quality was assessed using the Agilent RNA 6000 Pico chips (Agilent; PN- 5067-1513), ensuring a minimum RIN number of 7. Second, 10-μm sections were cut on the Visium Spatial Tissue Optimization Slide (PN-1000193) to assess optimal tissue permeabilization time. FInally, a 10-μm section was mounted on a Visium Spatial Gene Expression Slide and then stained for H&E staining and imaged using the NanoZoomer S210 (Hamamatsu; NP-C13239) to assess tissue morphology and quality. Following protocol instructions, the sections were then permeabilized for 18 min, then tissue was lysed, and reverse transcription was performed followed by second strand synthesis and cDNA denaturation. Spatially barcoded, full-length cDNAs were amplified by PCR for 16 and 17 cycles, depending on the initial concentration previously determined by qPCR. Indexed sequencing libraries were generated via end repair, A-tailing, adaptor ligation, and sample index PCR. Size distribution and concentration of full-length GEX libraries were verified on an Agilent Bioanalyzer High Sensitivity chip (Agilent). Finally, sequencing of GEX libraries was carried out on a NovaSeq 6000 sequencer (Illumina) aiming at approximately 40,000 pair-end reads per spot.

Data pre-processing

For the analysis of the spatial transcriptomic data, SpaceRanger software (10x Genomics, v.1.3.0) was used to map the sequenced reads, correct amplification bias and obtain the count matrix. The mouse genome (mm10) was used as the reference. The filtered feature expression matrices generated were then used as input for downstream analysis with Seurat⁷⁵ (v.4.4.0) in R (v.4.3.1).

Quality control and data normalization

To ensure quality of the data, spots not overlapping tissue were removed previously to the SpaceRanger mapping with the Loupe Browser software (10X Genomics). Quality metrics were calculated on a per-slide basis to preserve biological variability. Differences in total UMI across spots were adjusted by log-normalization using the NormalizeData() function from Seurat. This function divides the raw gene counts for each cell by the total counts of that cell and multiplies it by the scale factor (10,000), which is then log-normalized as log(1+x). Genes not expressed in any spot overlapping tissue were also removed.

Feature selection and dimensionality reduction

To annotate the distinct lung regions of each Visium slide, we used the FindVariableFeatures() function to extract the top 3,000 HVGs and capture major axes of biological variability. Data were then scaled with ScaleData() to z-score. Principal component analysis was performed, and the top 50 principal components were retained for subsequent analysis steps.

Clustering and annotation

To perform clustering, the FindNeighbors() function was applied together with the Leiden⁷⁶ community detection algorithm. Sample-specific resolutions were chosen, ranging from 0.2 to 0.6 in the function FindClusters(). Lastly, DEGs for each cluster were identified with FindAllMarkers() function and Wilcoxon rank-sum test. When needed, clusters showing high heterogeneity were sub clustered and markers re-calculated. Each region of the lungs was then annotated considering the DEGs together with the haematoxylin and eosin staining.

Downstream analysis

To estimate cell-type composition in each spatial transcriptomic spot, we performed deconvolution using a single-cell reference dataset from the LungMap project⁶², using the seeded non-negative matrix factorization (NMF) regression approach implemented in SPOTlight (v.1.0.3)⁷⁷. Spots with a predicted composition of neutrophils of ≥10% were annotated as neutrophil-enriched.

Signature scoring was performed using decoupleR⁷⁸ (v.2.8.0). In brief, the univariate linear model (ulm) approach was applied to compute similarity (enrichment) scores by testing the association between gene expression and the neutrophil hub signatures derived from our single-cell RNA-seq data, thereby quantifying signature activity within each Visium spatial transcriptomics spot.

To map neutrophil signatures from the spatial transcriptomic data onto NeuMap, we first carried out differential expression analysis across healthy, flu-infected, and cancer samples using the FindAllMarkers() function in Seurat (v.4.3), with a logistic regression framework. Genes were included if expressed in at least 25% of cells in one group and showed a minimum log-fold change of 0.25. Significance was defined by adjusted P values (Benjamini–Hochberg correction) with a threshold of P ≤ 0.05.

The top 50 significantly DEGs (ranked by log₂ fold change) were selected for module scoring using Seurat’s AddModuleScore() function. For each gene set, a module score was computed by averaging expression levels and comparing against a background of control genes with matched expression, thereby controlling for overall expression bias.

Analyses of PDAC and myocardial infarction (MI) models

Spatial transcriptomic data were analysed for the PDAC and MI mouse datasets^7,23 and used the Seurat package (v.5.1.0) in R, with three biological replicates included for each condition. Raw count matrices were first filtered to remove unexpressed genes. Low-quality spots were excluded on the basis of thresholds for the number of detected genes, total UMI counts, and mitochondrial gene content. Spots with abnormally low or high total counts, low gene detection, or mitochondrial percentages exceeding dataset-specific thresholds were considered low quality and discarded. We normalized using SCTransform^79,80. Highly variable features were identified using the FindVariableFeatures() function, selecting the top 10% of genes by variability within each dataset. Data were then scaled using ScaleData() to centre gene expression values and apply z-score transformation. Dimensionality reduction was performed using PCA via RunPCA() on the selected HVGs. Neighbourhood graphs were constructed with FindNeighbors() on the basis of the first 20 principal components, followed by clustering with the Leiden algorithm and a resolution parameter set to 0.5 using FindClusters(). Low-quality clusters lacking underlying tissue structure were identified and removed. UMAP embeddings were computed on the same 20 principal components using RunUMAP() for visualization.

For single-cell referencing, we used a publicly available dataset (GSE141017and GSE176092)^23,80. Quality control was applied by removing unexpressed genes and low-quality cells on the basis of gene counts, UMI counts, and mitochondrial gene content. The dataset was normalized using SCTransform, followed by identification of HVGs (top 10%), scaling, and PCA using 30 (MI) to 40 (PDAC) components. Batch correction across samples was performed using FindIntegrationAnchors() and IntegrateData(). Clustering was performed with the Leiden algorithm (resolution = 0.5) after computing neighbours using the top 30–40 principal components. UMAP was used for visualization. We annotated cell-type identities using SingleR (v.2.6.0) with reference profiles from the MouseRNAseq ImmGen databases via the celldex package (v.1.14.0) and curated marker genes from the dataset’s own clustering results. Marker genes for each cell type were extracted using FindAllMarkers(). Cell-type-specific markers were used in subsequent annotation and deconvolution steps. Cell-type deconvolution of spatial transcriptomic data were conducted using the SPOTlight package (v.1.0.3), using a seeded NMF regression approach. Spots with neutrophil compositions of 10% or higher were labelled accordingly.

For neutrophil hub signature scoring and subtype analysis, we identified the top 15 marker genes for each hub, and module scores were calculated using AddModuleScore(). Seeded K-means clustering (K = 7) was performed on neutrophil-labelled spots using the hub marker signatures. Clusters were annotated on the basis of the most specific and abundant signature. Clusters with no dominant signature were left unclassified. This procedure was repeated for macrophage and T cell subtype signatures^81,82. For MI-associated fibroblasts⁸³, a different strategy was used by classifying the cell subtypes in the reference single-cell dataset and predicting their respective abundances directly via deconvolution. Spatial annotations were derived by integrating information from the Leiden clusters, cell-type deconvolution, histological inspection, and expression of cancer-specific signatures. GSEA and over-representation analyses were used to characterize and differentiate distinct tumour cores. Finally, to investigate spatial relationships between spots, we constructed a graph using the igraph::graph() function (v.2.1.4) on the basis of a distance matrix computed from spot coordinates (stats::dist()). The graph was tuned to include only immediate neighbours on the basis of the 2D spatial grid structure.

Spatial analysis of human lung specimens using Visium HD

Human tissue microarray samples were used under protocol 2019-5253, which was reviewed and approved by the McGill University Health Centre (MUHC) Research Ethics Board, specifically by the MUHC co-Chair of the Comité d’éthique de la recherche du CTGQ panel. Human lung tissue specimens were obtained through protocols approved by the McGill University Health Centre Institutional Review Board (IRB 2014-1119). From these samples, tissue microarrays were constructed by a pathologist on the basis of intratumoural neutrophil abundance, using 1-mm cores sampled from FFPE pulmonary invasive adenocarcinomas with high-grade predominant solid architecture and adjacent non-tumorous lung tissue. Samples were derived from eight patients. Sections (5 μm) were mounted onto Visium CytAssist HD slides (10x Genomics) and processed following the manufacturer’s protocol. In brief, FFPE tissue sections underwent deparaffinization, decrosslinking, probe hybridization, ligation and extension, followed by spatial barcoding and sample indexing. Final library quality and fragment size were assessed using an Agilent Bioanalyzer High Sensitivity chip. Libraries were sequenced at the McGill University Genome Centre on an Illumina NovaSeq X Plus platform (1.5B reads, PE100 per lane). Spatial gene expression data were processed using Space Ranger (10x Genomics), and high-resolution spatial transcriptomic profiles were generated for downstream analysis.

FASTQ files from the Visium HD experiment were mapped to the GRCh38-2020-A reference genome using the Visium_Human_Transcriptome_Probe_Set_v2.0_GRCh38-2020-A probe set with Space Ranger (v.3.1.1, 10x Genomics). Sample areas were manually selected using Loupe Browser, and the registration file H1-RTF6MBB-A1-fiducials-image-registration.json was provided to spaceranger count via the –loupe-alignment argument. High-resolution H&E and CytAssist images were passed with the –image and –cytaimage arguments, respectively.

Cell-level transcriptomic profiles were reconstructed using the Bin2cell package⁸⁴ (v.0.3.2). First, an AnnData object was created with b2c.read_visium, then filtered to retain bins with at least one count and genes expressed in at least three spots (min_counts = 1, min_cells = 3). The H&E image was scaled using b2c.scaled_he_image with mpp = 0.38. Due to variability in bin sizing across the array, 2 μm bins exhibit slight differences in width/height, leading to a striped appearance. To correct this artefact, the b2c.destripe function was applied. Nuclei segmentation was then performed on the H&E using b2c.stardist with parameters: stardist_model = “2D_versatile_he”, prob_thresh = 0.1, and nms_thresh = 0.1. Segmented nuclei were expanded using b2c.expand_labels with algorithm = “volume_ratio”. To recover additional cells not segmented via H&E, the Stardist fluorescence model was applied to a σ-smoothed gene expression image generated using b2c.grid_image (mpp = 0.38, sigma = 5). b2c.stardist was re-run with stardist_model = “2D_versatile_fluo” and the same thresholds. Cells not segmented with the H&E-based model were assigned secondary labels from the fluorescence model. Finally, bins were grouped into cells with the b2c.bin_to_cell function. To reconstruct the Segmentation Polygon Mask we converted the.npz mask output from Bin2cell into a data frame using pd.DataFrame.from_dict (pandas v.2.2.3). This file was then processed in R (v.4.4.1) using the concaveman package (v.1.1.0) to generate cell polygons for visualization.

SingleCellExperiment⁸⁵ objects were created in R by importing the Bin2cell.h5ad files with the h5ad2sce function from the schard package (v.0.0.1). Low-quality cells were removed if they met any of the following criteria: fewer than 10 counts, fewer than 10 unique genes, area <8 µm², or >15% mitochondrial reads. Additionally, cells labelled only through gene expression segmentation (secondary labels) were excluded upon manual inspection for falling outside tissue boundaries. In a second quality control step, the isOutlier function from the scuttle package (v.1.16.0) was used to flag outliers in transcript count and transcript density (log₂(counts/area)) using parameters nmads = 2.5, type = “both”, and log = FALSE. Identified outliers were excluded from downstream analysis. For pre-processing, quality-filtered SingleCellExperiment objects were normalized using logNormCounts from scuttle. HVGs were identified with modelGeneVar (scran⁶ v.1.34.0), blocking by patient ID (block = sce$patient). HVGs were selected using getTopHVGs (fdr.threshold = 0.05). PCA was performed using runPCA (scater v.1.34.0)⁸⁶ on the HVGs (subset_row = hvgs), and the first 15 principal components were retained on the basis of elbow plot inspection. UMAP dimensionality reduction was computed using runUMAP (scater).

To cluster cells, we built a shared nearest-neighbour graph using buildSNNGraph from scran with type = “jaccard” and use.dimred = “PCA”. Louvain community detection was performed using cluster_louvain (igraph⁷ v.2.1.1) at multiple resolutions. DEGs between clusters were identified using findMarkers (scran) with direction = “up”. Top-ranked markers for each cluster were selected from the ‘Top’ column of the results and inspected manually to guide spatial annotation. Hallmark gene sets (for example, angiogenesis, hypoxia) were obtained from MSigDB using the msigdbr package (v.7.5.1) with species = “Homo sapiens” and category = “H”. Immune cell-type signatures⁸⁷ were retrieved using category = “C8” and filtered to retain adult lung signatures while excluding fetal profiles. Signature scoring was performed using AUCell (v.1.28.0). Gene expression rankings per cell were computed with AUCell_buildRankings, followed by AUC calculation with AUCell_calcAUC. For visualization, cluster-level mean AUC scores were obtained using aggregateAcrossCells (scuttle⁸⁶) with statistics = “mean” and use.assay.type = “AUC”.

To assign neutrophil clusters to NeuMap hubs, human gene signatures corresponding to spatial clusters C1–C5 were first converted to mouse orthologs using homology data from the Mouse Genome Informatics (MGI) database (https://informatics.jax.org). We then computed gene module scores for each hub using the AddModuleScore() function in Seurat. Module scores were calculated by averaging the expression of hub-associated genes and comparing them to control gene sets with similar expression levels, thereby normalizing for baseline expression and minimizing bias. For visualization, enrichment scores were averaged by hub across clusters and scaled per signature to enable comparison in the resulting heat map.

Neighbourhood analysis

Spatial transcriptomic data were analysed using Seurat (v.4.3) in R (v.4.2). A spatial proximity graph was constructed by computing the k-nearest neighbours (k = 6) from each spot’s x,y coordinates using the RANN package (v.2.6.2). Edges exceeding 250 units in Euclidean distance were excluded to account for realistic cell–cell interaction radii. An undirected graph was generated using the igraph package, with edge weights corresponding to physical distances. We implemented two complementary functions to quantify spatial cell-type context: (1) neighbourhood frequency analysis. For a given set of target spots (for example, neutrophils), their first-order neighbours were identified within the spatial graph. For each neighbour, the cell-type label was extracted, the number of neighbouring cells of each type was counted per target spot and aggregated across all targets to calculate the mean and 95% confidence interval for each cell type. This allowed identification of the most frequently co-localized cell types around a given population. (2) neighbourhood composition by Hub: In a separate analysis, both the target spot and its neighbours were pooled to represent a local “neighbourhood.” Cell types within each neighbourhood were classified using a curated cell-type annotation. For each neutrophil hub, cell-type counts were summed and normalized to percentages, excluding spots annotated as ‘unassigned’ to avoid skewing proportion estimates. This enabled comparative analysis of cellular composition across microenvironments.

MACSima imaging cyclic staining

Sample preparation and image acquisition

Multiplex immunohistochemistry of lungs from naive, flu-infected, or tumour-induced mice was performed using a MACSima imaging system (Miltenyi Biotec). In brief, cyclic immunofluorescence imaging consisting of repetitive cycles of immunofluorescent staining, sample washing, multi-field imaging, and signal erasure by photobleaching was performed. Cryosectioned fixated lungs from the 3 groups were placed on microscopy slides and MACSwell sample carriers were mounted and blocked using a blocking buffer containing 10% BSA and 2% goat serum for 1 h at room temperature before lungs were preincubated with an antibody to IFIT1 (ab236256, Abcam, 1:100) overnight at 4 °C. Thereafter, nuclei were counterstained with DAPI before samples were placed in the MACSima imaging system. Neutrophil subsets were identified using the following antibodies: anti-CD11b-APC (clone M1/70.15.11.5, Miltenyi Biotec, 130-113-239; 1:50); anti-CD45-FITC (clone REA737, Miltenyi Biotec, 130-110-658; 1:50); anti-CXCR2-PE (clone SA044G4, BioLegend, 149303; 1:50); anti-Ly6C-PE (clone REA796, Miltenyi Biotec, 130-111-916; 1:50); anti–MHC-II-APC (clone REA813, Miltenyi Biotec, 130-112-388; 1:50); anti-PD-L1-APC (clone 10 F.9G2, BioLegend, 124312; 1:50); anti-CD14-PE (clone Sa14-2, BioLegend, 150106; 1:50) Anti-IFIT1 (Polyclonal, Abcam, Ab70023; 1:50). Additionally, a conjugated anti-rabbit antibody (polyclonal, Sigma-Aldrich, F9887; 1:100) in the first cycle to identify IFIT1. The tissue location was characterized using anti-Podoplanin-PE (clone 8.1.1, BioLegend, 127408; 1:50).

Data analysis and visualization

Images were stitched and pre-processed using MACS iQ View Analysis Software (Miltenyi Biotec, v.1.3.2) and representative overlay pictures were displayed. For downstream analysis, cells were segmented on the basis of the DAPI signal using the StarDist plugin in ImageJ (US National Institutes of Health) and the donut algorithm in MACS iQ View. The segmented data were then exported for further analysis to FlowJo (BD Biosciences), where neutrophils were identified as LY6G^hi cells using a threshold value. In flu-infected lungs, neutrophils were characterized as intravascular or extravasated on the basis of the podoplanin signal in the segmented neutrophils. Neutrophils from all ROIs were concatenated and phenotypically analysed by dimensional reduction using the UMAP plugin and unsupervised clustering using the FlowSOM plugin. FlowJo and R were used for visualization.

Mathematical modelling for blood neutrophil diagnosis

To model the predictive value of the distribution of blood neutrophils in NeuMap, we used the density overlap represented by the Bhattacharyya index. We favoured density versus spatial overlap to reduce the impact of outliers in our calculations. We also favoured using ten selected regions of NeuMap over the seven hubs delineated in Fig. 1e to gain spatial resolution of cell distribution over the specific areas of NeuMap where blood neutrophils from the tested conditions tended to concentrate. In separate analyses, we found that the spatial overlap was not based on densities, and the use of seven hubs provided significantly less resolution in the overlap barcodes (not shown).

Density state estimation

The UMAP coordinates \(\x_j,y_j\\) of the neutrophils collected from the tissue hubs and blood samples were pipelined into a computational approach to estimate the probability density functions (PDFs) over the Neumap transcriptomic space associated with these datasets. The kernel density estimators \(\hatf\,(x,y)\) were of the form⁸⁸,

where K is the Epanechnikov kernel

with n, h_x and h_y denoting the number of points and the bandwidths along two orthogonal directions, respectively. This kernel choice, which is symmetric and normalized along each orthogonal direction, ensures a smooth, bounded and efficient computation of the density estimates. The determination of each kernel’s bandwidths, which influences the smoothness and accuracy of the resulting PDF estimate, was made via the Sheather–Jones method⁶⁵. This technique uses a data-driven approach that minimizes the mean integrated square error of the estimated density function. By iteratively adjusting the bandwidths h_x and h_y and evaluating their performance, Sheather–Jones effectively balances bias and variance, resulting in accurate density estimates, particularly for tissue hubs and blood samples that present non-Gaussian distributions or multimodality.

Overlap integration

Upon estimating the PDFs of our datasets, we assessed the degree of overlap between pairs of these functions. For this purpose, we used the Bhattacharyya index, which yields a measure of the amount of overlap between two PDFs \(f(x,y)\) and \(g(x,y)\). The Bhattacharyya index was defined by the integral expression

where the double integral runs over the entire domain defined by NeuMap. The Bhattacharyya index gives a scalar value in the interval [0, 1], enabling a direct interpretation of the overlap: values near 1 suggest a high degree of similarity, or almost perfect overlap, between the two distributions, indicating that the neutrophil states are nearly indistinguishable within the dimensionally reduced transcriptomic space. Conversely, values near 0 denote little to no overlap, pointing to distinct neutrophil states with minimal similarity in their respective distributions within the dimensionally reduced transcriptomic space.

To calculate the Bhattacharyya index, we used an adaptive Monte Carlo method⁸⁹ that combines adaptive importance and stratified samplings over multiple iterations, thus optimizing the sample distribution around the peaks of the PDFs and thereby reducing the standard deviation in the estimates. This approach yielded accurate values of the multidimensional integrals and hence offered robust measures of the overlap between the PDFs. While the Bhattacharyya index does not constitute a true probability measure, its bounded nature makes it a valuable score for comparing the similarity of data distributions in a normalized manner. From the computed Bhattacharyya indexes, the resulting barcodes were generated. To perform our calculations, we used the R programming language to analyse all datasets. Specifically, we used the MASS, graphics, stats and vegas R packages.

Estimation of neutrophil lifetimes

To quantify the neutrophil half-lives and transit times, we used an age-structured mathematical model as previously proposed¹. This model effectively captures the temporal variation in the proportion of labelled neutrophils following the administration of the BrdU pulse. Let \(u=u(t,a)\) denote the density of neutrophils, which at time t, have an age a. We assume that their age ranges in the interval \(a\in [0,a_\max ]\), where \(a_\max \) is the maximum age (or the maximum lifespan) a neutrophil can achieve in the different examined tissues. In practice, this age can be taken sufficiently large without appreciably altering the numerical results for the entire neutrophil population. To describe the temporal dynamics of the age distribution in neutrophils, we considered the following first-order linear transport partial differential equation

The left-hand side of equation (1) represents the temporal change in the number of neutrophils along with their corresponding age. The first term on the right-hand side accounts for neutrophil death. The death time, \(\tau (a)\), generally depends on the age of the neutrophil. The introduction of a flux function, \(\phi (t,a)\), encapsulates the net recruitment of neutrophils entering or leaving the target tissue. Assuming that at time \(t=0\) no BrdU⁺-labelled neutrophils of any age have yet arrived at tissue i, the initial condition is \(u(0,a)=0\). Therefore, the exact solution to (1), obtained using the method of characteristics for first-order partial differential equations, is

Equipped with equation (2), we computed the total number of neutrophils at time t and tissue i, irrespective of their age, via the following integral

To connect equation (2) with the different scenarios addressed in the experiments, the net flux \(\phi _i(t,a)\) corresponded to one synchronous wave of neutrophils after administration of the BrdU labelling. The chosen functional forms for \(\tau _i(a)\) and \(\phi _i(t,a)\) in our modelling were

where τ_i, α_i, t_i, and σ_i are constant parameters. Inserting equations (2) and (4) into equation (3), we arrive at the following exact formula for the total number of neutrophils at a given tissue i at time t

where

is a normalization parameter for tissue i, and \(\mathrmErf(x)=\frac2\sqrt\rm\pi \int _^x\rme^-\xi ^2\rmd\xi \) is the error function. The parameters τ_i, m_i, t_i, and σ_i for blood, bone marrow and spleen were computed via a nonlinear regression analysis using the corresponding time series measured from day 1 to day 7 for each tissue compartment. Once these four parameters were found for each tissue i, the mean half-life time \(t_1/2^(i)\) and transit time \(t_\mathrmtran^(i)\) were estimated from the normalized profile (equation (5)). To do that, the transit time was identified as the time at which this unimodal profile achieves its maximum (100% of the BrdU-labelled neutrophils in tissue i). Subsequently, the 50% level was set as a reference for the mean half-life time \(t_1/2^(i)\). The approach for calculating these two lifetimes is illustrated in Extended Data Fig. 7f. The bands shown in this figure correspond to a confidence level of 0.75.

To carry out the density state estimation and overlap integration, we used the R programming language to analyse all datasets. Specifically, we used the MASS v.7.3.61, graphics v.4.4.3, stats v.4.0.3 and v.4.4.3., and vegas v.2.1.4 R packages. The nonlinear regression and statistical analysis were performed with Matlab (R2024a) using the functions fitnlm and coefCI.

Quantification and statistical analysis

Data from experiments are represented as mean values ± s.e.m. All parameters analysed followed normal distribution as tested by D’Agostino–Pearson test unless indicated in the figure legend. Unpaired two-tailed t-test was used when two groups were compared, and comparison of more than two datasets was done using one-way analysis of variance (ANOVA) with Tukey’s post-test or two-way ANOVA. Log-rank analysis was used for Kaplan–Meier survival curves. Statistical analysis was performed using GraphPad software. Statistics on the RNA sequencing are indicated in the analysis section. A P value below 0.05 (*) was considered statistically significant; P ≤ 0.01 (**) and P ≤ 0.001 (***), as well as nonsignificant differences (NS), are indicated accordingly.

Reporting summary

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