Study design
The overarching goal of this study was to characterize innate memory of NAMs in heterologous infections as well as to investigate NAM intrinsic mechanisms that drive disease tolerance. To do so, we used a combination of genomic approaches, in situ microscopy, flow cytometry as well as functional assays in an in vivo murine system.
Mice
All mice used in this study were bred and maintained in the Alexandria Center for the Life Sciences animal vivarium or purchased from commercial vendors when mentioned, and were males aged between 7 and 16 weeks at the start of the experiment. We used 3 main genotypes throughout the study: WT mice C57BL/6 mice were either purchased from Charles River Laboratories (C57BL/6NCr-NIH, 556) and left to be acclimatized in the NYU vivarium for 2 weeks or were bred in-house in the NYU Vivarium; (2) NAM-DTR (Cd169-cre mice crossed to B6N.129P2-Cx3cr1tm3(DTR)Litt/J mice) (as described in the supplementary information of ref. 3) were maintained on a C57BL/6 background and (3) Arg1-cKO mice (CD169-cre mice crossed to C57BL/6-Arg1tm1Pmu/J mice) were purchased from Jackson Laboratories. An F1 cross generated CD169-cre × Arg1fl (homo/homo), further crossed to Arg1fl to generate an experimental CD169-cre × Arg1fl (het/homo) was developed in house and maintained on the B6 background and are referred to as Arg1-cKO mice in this study. The numbers of mice per experimental group are indicated in the figure legends. Throughout the study, all mice were maintained with food and water ad libitum under a 12 h dark–light cycle, with 50–60% humidity and at 20–25 °C in a pathogen-free facility at New York University Langone Health Center. All of the experiments were performed with approval by the New York University Langone Health Center Institutional Animal Care and Use Committee and in accordance with guidelines from the National Institutes of Health, the Animal Welfare Act, and the US Federal Law.
Parasite culture and inoculation of mice
N. brasiliensis L3 larvae were cultured and maintained as previously described43 in a Petri dish culture containing 50% charcoal and 50% spaghnum peat moss. The larvae were isolated from cultures using a modified Baermann apparatus with 400 U penicillin, 400 µg ml−1 each of Neomycin (Gibco) and streptomycin in sterile 1× PBS, then washed with sterile PBS three times. Mouse colonies were maintained by culturing faecal and caecal palettes from either Stat6−/− or C57BL/6 mice. To maintain similar worm infectivity, L3-stage juvenile larvae were cultured from plates within 8 weeks of preparation. For all of the experimental infections, mice were subcutaneously inoculated with 200 µl suspension of 550–600 N. brasiliensis L3 by lifting the skin on the head between both ears under isoflurane anaesthesia. For experimental consistency of infection based on worm motility, mice were infected with L3-larvae isolated from plates that had been prepared within 8 weeks to maintain the same degree of infectivity.
Influenza infection
Acute respiratory viral infection was introduced by an i.n. administration of influenza virus (H1N1/PR8), which was purchased from Charles River Laboratories. Mice were infected with 75 to 150 EID50 in a 30 µl suspension with sterile PBS through an i.n. route under anaesthesia and analysed at various timepoints as indicated in the main text.
Lung macroscopic analysis and histopathology. After euthanasia with CO2, mouse lungs were collected, fixed in 4% paraformaldehyde (PFA) for 72 h at 4 °C on a shaker. Fixed lungs were washed with PBS three times, then transferred into 70% ethanol, and paraffin-embedded for sectioning. Embedded lungs were sectioned at 5 µm and stained with H&E. Lung damage was scored on five randomly selected high-power fields (×40) based on neutrophil infiltration and aggregation into the lung parenchyma using previously described methods29,57.
Nucleus detection and quantification using QuPath
Image analysis for histochemical stains was conducted on QuPath (v.0.6.0-x64, macOS) to quantify nuclear profiles. ×40 histology images were scanned and uploaded at the NYU core on Omero. High-resolution images were imported into QuPath project and pixel calibration (μm per pixel) was confirmed from metadata and/or a scale bar to enable area-based normalization. Tissue regions of interest (ROIs) were defined using annotations to include tissue and exclude background and artifacts (folds/tears). For H&E slides, stain vectors were estimated, and nuclei were detected using QuPath’s watershed-based cell detection with the haematoxylin channel as the detection image; key parameters (threshold and min/max nucleus size constraints) were optimized on representative images and then applied uniformly across all images/conditions. Nucleus counts were computed as the number of detections within each ROI, and the nucleus density was reported as the total numbers within the ROI. Detection overlays were visually inspected for quality control, and ROIs with major artifacts were refined or excluded before exporting final counts and densities.
Gross pathology
Before collection, lung images were collected using a Samsung Galaxy S22+ smartphone and high-resolution images were uploaded to a computer, then exported to illustrator for annotation and to highlight disease and morphological features of PR8-infected lungs.
Tissue preparation for immunofluorescence and confocal microscopy
Mice were euthanized using CO2 exposure and right lung lobes were collected and subsequently fixed as previously described3. Tissues were fixed in PFA, l-lysine, sodium periodate buffer (PLP, 0.05 M phosphate buffer, 0.1 M l-lysine, pH 7.4, 2 mg ml−1 NaIO4 and 40 mg ml−1 PFA) overnight at 4 °C on a shaker. The next day, tissues were dehydrated in 30% sucrose overnight at 4 °C and subsequently embedded in optimal cutting temperature compound (OCT) medium. Frozen tissue sections were sectioned using the Leica CM1850/Leica 3050S system at a thickness of 20 μm. Fc receptors were blocked with anti-CD16/32 Fc block antibody (BioLegend) diluted in PBS containing 2% serum and 2% FBS for 1 h at room temperature. Fc receptor was blocked with anti-CD16/32 Fc block antibody (93, BioLegend) diluted in 1× PBS containing 2% donkey or goat serum, 2% FBS and 0.01% Triton X-100 for 1 h at room temperature. The sections were stained with Ly6G-BV421 or Ly6G-PE or Ly6G-BV510 (1A8, BioLegend), CD169-ef660 (Ser-4, Invitrogen), anti-histone H3 (Abcam), Siglec-F-BV421 or Siglec-F AF488 or Siglec-F PE (E50-2440, BD Biosciences), IL-5-BV421 (TRFK5, BioLegend), keratin 5 (polyclonal, BioLegend), SFTPC (polyclonal, Sigma-Aldrich) and Acetylated α-tubulin (6-11b-1, ThermoScientific) and diluted in blocking buffer for 1 h at room temperature. The sections were then washed with 1× PBS and subsequently stained with donkey anti-goat 488 (Thermo Fisher Scientific), goat anti-rabbit 488 (Thermo Fisher Scientific), goat anti-rabbit 546 (Thermo Fisher Scientific), goat anti-mouse 488 (Thermo Fisher Scientific) in blocking buffer for 1 h at room temperature and then washed with 1× PBS. After the final wash, the slides with sections were mounted using Immun-mount mounting medium (Thermo Fisher Scientific) and cover glasses with a thickness of 0.13 to 0.17 mm (Thermo Fisher Scientific) were placed and then left to dry overnight in drawers; the slides were stored in the dark at 4 °C until ready for microscopy analysis. Fluorescence was detected using the Zeiss LSM 880 confocal microscope (Carl Zeiss) equipped with 405, 488, 514, 561, 594 and 633 nm solid-state laser lines, a 32-channel spectral detector (409 to 695 nm) and ×10/0.3, ×20 Plan-Apochromat 0.8, ×40 and ×63/1.40 objectives. The Zen Black (Carl Zeiss) software suite was used for data collection. Whole-section imaging was performed using an All-in-One Fluorescence Microscope BZ-X800 (Keyence) equipped with an eGFP, Cy3/TRITC and Cy5 filter cube and a ×20 Plan Fluor objective. Keyence images were stitched using the BZ-X800 analyser. Confocal imaging data were processed and analysed using Imaris software v.8.3.1 (Bitplane; Oxford Instruments).
Tissue homogenization
Lungs were collected at the designated timepoints and washed in 1× PBS before being placed into 2 ml tubes containing 500 µl of sterile 1× PBS and one 2 mm ceramic bead (MP Biomedical). Lung tissue was homogenized using the FastPrep24 instrument at 4.0 m s−1 for 40 s. After homogenization, lung debris was centrifuged at 15,000g for 5 min. The supernatant was aliquoted in 100 μl aliquots and frozen at −80 °C until being used for analyses.
Preparation of cell suspensions for flow cytometry and cell sorting
Single-cell suspensions of lung tissues were prepared for flow cytometry. Lung tissues were injected and incubated with RPMI 1640 (Lonza) medium containing Liberase (5 mg ml−1; Sigma-Aldrich for myeloid cell evaluation) or collagenase I (100 U ml−1 for lymphocyte evaluation) enzyme supplemented with 10% FBS (Gibco), 0.2% CaCl2, 0.2% MgCl2, deoxyribonuclease I and 1% HGPG (1 mM HEPES, 5 mM l-glutamine, penicillin–streptomycin (10,000 U ml−1) and gentamicin (5 μg ml−1) (pH 7.5)) for 30 min at 37 °C. Enzyme digestion was inactivated by addition of PBS containing 1 mM EDTA and 2% FBS. Lung tissue was dissociated into single-cell suspensions, and red blood lysis was performed using RBC lysis buffer (Alfa Aesar). Lung cells were resuspended in fluorescence-activated cell sorting (FACS) buffer (1× PBS (Gibco), 2% FBS, and 1 mM EDTA). Fc receptors were blocked with anti-CD16/32 Fc block antibody (BioLegend) and stained with the indicated antibodies (below) for 30 min at 4 °C. Cells were fixed with 2% PFA for 10 min at room temperature and resuspended in FACS buffer. For intravascular staining to label blood versus tissue-resident immune cells, 1 µg of anti-CD45 BUV661 (30-F11, BioLegend) diluted in 150 µl of PBS was injected retro-orbitally into mice under anaesthesia. Mice were euthanized 3–4 min later, and organs were collected for analyses. The cell suspension was processed on the Bio-Rad ZE5 instrument (Bio-Rad), and data were analysed using FlowJo software (BD Biosciences). For adoptive transfer experiments, 3 mice each from the PBS-treated or N. brasiliensis-infected mice were pooled after homogenization and before staining. Cells were stained with DAPI for continuous labelling of dead or dying cells in addition to staining for CD45, Lin− dump channel (CD3, NK1.1, B220, CD31, EPCAM), CD11b, CD11c, MerTK, CD64 and Siglec-F. Cells were sorted using the Symphony S6 SORP Sorter at NYU Flow Cytometry Core Facility. The post-sort purity was checked for each sample and maintained at 95–98%. Aseptic conditions were maintained and sorted NAMs were washed twice to remove FBS using a swinging-bucket centrifuge at 4 °C. Washed cells were reconstituted in sterile 1× PBS at 1 × 104 cells per 80 μl for intratracheal administration.
Antibody stains. Fc receptors were blocked with anti-CD16/32 Fc block antibody (93, BioLegend) and evaluation for lung myeloid cells was done by staining with Live Dead UV (AF350 NHS Ester, Thermo Fisher Scientific) and staining for CD45 (MEL-14, BioLegend). For the myeloid-specific panel, we used the following antibodies: Siglec-F (E50-2440, BD Bioscience), MHC-II (M5/114.15.2, BioLegend), CX3CR1 (SA011F11, BioLegend), F4/80 (BM8, BioLegend), CD11b (M1/70, BD Bioscience), Ly6C (HK1.4, BD Bioscience), MerTK (DS5MMER, Thermo Fisher Scientific), Ly6G (IA8, BioLegend), CD11c (N418, Thermo Fisher Scientific), CD64 (X54-5/7.1, BioLegend), CD169 (Ser-4, Thermo Fisher Scientific), CCR3 (J073E5, BioLegend). For the lymphocyte panel, we used the following antibodies: CD19 (6D5, BioLegend), B220 (RA3-6B2, BioLegend) CD4 (RM4-5, BioLegend), CD3 (17A2, Thermo Fisher Scientific), CD44 (IM7, BioLegend), CD11a (M17/4, BioLegend), CD62L (MEL-14, BioLegend), delta TCR (GL3, BD Biosciences), CD127 (A7R34, BioLegend), NK1.1 (PK136, BioLegend); CD8a (53-6.7, BioLegend) to characterize lymphocyte populations. Additionally, we used CD8+ specific H2-Db, influenza A NP 366–374, ASNENMETM (PE-labelled) or NP-Tet 311–325 (APC-labelled) for analysis of influenza-specific CD8+ T cells at 6 days after infection.
Stained cells were analysed on the Bio-Rad ZE5 instrument with Everest Software, and data were analysed annotated using FlowJo software. For cell sorting, lung tissue was treated in the same way without fixation. After staining, cell pellet was directly resuspended in 1× PBS with 1 mM EDTA and 2% FBS. Cells were sorted on the BD FACS Aria II instrument or BD Symphony S6 SORP Sorter.
In vivo cell depletions
In vivo depletion of NAMs was done locally by i.n. administration of 15 ng per mouse of DT (50 µl) twice 3 days apart or systemically by administration of a single i.p. injection of 40 ng g−1 (200 μl). Control mice were injected with 200 µg of anti-IgG2b (LTF-2, BioXcell) i.p. (200 µl). In vivo depletion of resident eosinophils (CCR3+Siglec-F+Ly6G−CD11b+Ly6C−) cells in vivo was done by injecting mice with 200 µg per mouse of anti-CCR3 (6S2-19-4, BioXCell) i.p. (200 μl). In vivo depletion of IL-4 was done using anti-IL4 neutralizing antibody (11B11, BioXcell) along isotype control anti IgG2a (RG7/1.30, BioXcell) and mice were treated with 200 µg per mouse (i.p.) once on day 7 after N. brasiliensis infection. In vivo depletion of IL-5 was done using anti-IL5 neutralizing antibody (TRFK5, BioXcell) along isotype control anti IgG1 (HRPN, BioXcell) and mice were treated with 200 µg per mouse (i.p.) once on day 7 after N. brasiliensis infection.
Recombinant protein treatments
Recombinant IL-4 was administered by formulating rIL-4 complex which is 1.5 µg per mouse of rIL-4 (11B11, Invitrogen, MA5-47800) plus 7.5 µg per mouse of anti-IL-4 (BioXcell, 6S2-19-4) administered i.p. (200 µl). Recombinant IL-5 was administered at a dose of 50 ng per mouse (R&D Systems, 405-ML-005/CF).
Adoptive transfers and challenge with PR8
For adoptive transfer experiments, 3 mice from each of the conditions of PBS-treated or N. brasiliensis-infected cohorts were collected and digested with Liberase as described earlier. Single-cell suspensions were stained for immune markers and subsequently sorted on the Symphony S6 Sort system. After a purity check, cells were resuspended in sterile 1× PBS after two washes, then reconstituted at a concentration of 1 × 104 cells per 80 μl. The mice were anaesthetized with ketamine + xylazine, then mounted onto slant plank and hung by their frontal teeth and administered 80 μl of cell suspension intratracheally. Mice were monitored for additional 30 min to ensure waking up from anaesthesia and no asphyxiation by cell transfer. Then, 14 h later, all of the mice were infected with 75 EID50 (sublethal dose) i.n. and monitored for weight loss daily until the end of study.
Cytokine and chemokine analysis
The supernatants from whole-lung lysates were used for determination of chemokines and cytokines using Bio-Plex Pro mouse cytokine 23-plex assay (Bio-Rad, M60009RDPD) or mouse IFN-alpha/IFN-beta kit (ProcartaPlex, EPX02A-22187-901) or viral titres.
Measurement of lung viral titres using plaque assay
Supernatants from whole-lung homogenates were collected at 1, 3 and 6 days after infection. Virus titres were determined on Madin–Darby canine kidney (MDCK) cells using a standard plaque assay. MDCK cells (1.5 × 106 cells per cm2) were seeded in six-well culture plates and incubated at 37 °C under 5% CO2 for 24 h. Cells were washed twice with 1× PBS before adding 200 μl of a serial tenfold dilution of mouse lung extracts (dilution prepared in 1× PBS and 0.3% BSA). After 1 h of incubation at 37 °C, the inoculates were removed from the cells and overlay medium (2 ml per well; Dulbecco’s modified Eagle’s medium/F-12 medium; 2 mM glutamine, 10 mM HEPES, 0.15% NaHCO3, 1% penicillin–streptomycin, 0.2% BSA, 0.6% Oxoid purified agar, 0.01% DEAE–dextran hydrochloride and TPCK-trypsin (1 μg ml−1)) was added to the cells, followed by incubation at 37 °C in 5% CO2 for 72 h. Subsequently, the cells were fixed with 4% PFA in 1× PBS (1 ml per well) for at least 4 h at room temperature. Plaque overlay was removed, and cells were washed with 1× PBS and set to dry overnight. Cells were then counterstained with 1% crystal violet for 1 min, washed briefly in water, and set to dry before plaque enumeration.
WBP analysis
Plethysmography tests were performed using previously established protocols58 with necessary modifications for this study. Naive, trained or influenza-infected mice were introduced to a plethysmography chamber (Buxco Systems), which was placed in an unlit biological safety cabinet, and the mice were acclimatized to the chamber for 10 min. After acclimatization, the mice were actively measured for a 10 min period. Mice were randomly assigned to the plethysmography chamber for the intended recordings. During the 10 min measurement period within the plethysmograph, the plethysmograph continually measured respiratory responses, with 2 s summaries for 11 various respiratory parameters such as f, Tdv, Mvb, Penh, Rpef, PIF, PEF, Ti, Te, EF50 and Tr. For each condition, we used 4–5 animals and computed the geometric mean ± s.e.m. For each parameter, calculations were made using Finepointe software (v.2.9.0) and exported to Excel for generating data graphs. Data in Figs. 1 and 4 are shown only for selected members chosen as described in other published studies with differences in f, Penh, Tp, EP50, Rpef.
RNA extraction and RNA-seq analysis
Sorted NAMs from lung cell suspensions (15–20,000 cells per cohort) from 12–15 control C57BL/6 mice or N. brasiliensis-trained mice were pooled for the RNA-seq analysis. Cells were obtained and sorted as indicated above. The RNAeasy plus micro kit (Qiagen) was used to extract total RNA from the sorted cells. RNA-seq libraries were prepared using the Illumina TruSeq Stranded mRNA Library Kit, and the individual libraries were pooled together. The samples were run over four lanes on the Illumina HiSeq2500 system with 101 bp paired-end reads.
Sequencing reads were mapped to the reference genome (mm10) using the STAR aligner59 (v.2.5.0c). Alignments were guided by a gene transfer format (GTF) file. The mean read insert sizes and their standard deviations were calculated using Picard tools (v.1.126) (http://broadinstitute.github.io/picard). The read count tables were generated using HTSeq60 (v.0.6.0), normalized based on their library size factors using DEseq261 and differential expression analysis was performed. The reads per million (RPM)-normalized BigWig files were generated using BEDTools62 (v.2.17.0) and bedGraphToBigWig tool (v.4). Gene set enrichment analysis was performed using the GSEA tool63. To compare the level of similarity among the samples and their replicates, we used two methods: PCA and Euclidean distance-based sample clustering. All of the downstream statistical analyses and plot generation were performed in R environment (v.3.1.1) (https://www.r-project.org/).
Differential expression analysis and PCA. The raw counts were provided as multiple different.txt files—one for each replicate. Therefore, we first merged them into a single counts matrix using the data.table (v.1.17.8) package. Next, the data were prepared for differential expression analysis using DESeq2 (v.1.48.2). Genes with low overall counts were filtered, samples were annotated based on training status and cell type, and a DEseqDataSet (DDS) object was constructed using the resulting count matrix and experimental design. To assess sample-to-sample variation and potential outliers before differential expression testing, raw count data in the DDS were transformed using DESeq2’s variance-stabilizing transformation, which reduces the mean-variance dependence inherent to count data. PCA was then performed using DESeq2’s plotPCA function and plotted using ggplot2 (v.4.0.1). Differential expression analysis was then performed using DESeq2. Gene-wise differential expression between trained and naive NAMs was assessed using the DESeq2 Wald test within a negative binomial linear model framework. P values were adjusted for multiple testing using the Benjamini–Hochberg method.
ATAC–seq sample prep and libraries
A total of 40,000–50,000 purified macrophages per replicate was washed once in cold sterile PBS and resuspended in 20–50 µl of cold lysis buffer (10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2 and 0.1% IGEPAL CA-630). Cell lysates were centrifuged at 500 g for 30 s to 1 min at 4 °C, and nuclei were resuspended in 25 µl of water. The transposition reaction mix consisted of 25 µl of TD buffer, 2.5 µl of Tn5 transposase (Illumina, 20034198) and 22.5 µl of nuclease-free water containing the nuclei. The mix was incubated for 30 min at 37 °C. Transposed DNA fragments were purified using the Qiagen Reaction MiniElute kit (Qiagen). Clean DNA was run with PCR amplification to incorporate the indexes, which were a universal index and a unique i7 index per sample (IDT primers were outlined previously64) and amplified by PCR for 11 cycles using the NEBNext high fidelity 2× PCR master mix (New England Biolabs). PCR products were cleaned with 0.9× volume of Ampure XP beads (Beckman Coulter), and the amplified fragment size was verified on the Agilent 2100 Bioanalyzer system (Agilent) with high-sensitivity DNA screentape. ATAC–seq libraries were quantified by qPCR using the Roche-kapa Illumina quantification kit. Paired-end sequencing was performed on the Illumina Novaseq 6000 system with a SP100 cycle kit sequenced in paired-end 50 bp format.
GSEA. Before performing pathway-level analysis all genes tested for differential expression were ranked by the DESeq2 Wald statistic. The MSigDB Hallmark gene sets for Mus Musculus were then obtained using the msigdbr package (v.25.1.1). Finally, the clusterProfiler package (v.4.16.0) was used to perform the GSEA analysis, where the pre-ranked list of genes tested for differential expression was analysed against the MSigDB mice hallmark pathways. Gene sets of size 10–500 were included, and the results were filtered using a FDR threshold of 0.25.
Motif enrichment analysis. To identify transcriptional regulators associated with the enriched pathways, transcription factor motif enrichment analysis was performed using HOMER (https://find.motifs.pl) (v5.1). The foreground set consisted of 543 leading-edge (core) genes derived from significantly enriched, upregulated MSigDB Hallmark gene sets identified by GSEA. As a background, 14,900 genes with cumulative expression of 30 TPM across 6 bulk RNA-seq samples were used to restrict the analysis to expressed genes. Motif enrichment was assessed in promoter regions spanning −300 bp upstream to +50 bp downstream of the annotated transcription start site, using default HOMER parameters and the mm10 mouse genome.
ATAC-seq pre-processing and analysisPreprocessing of reads was done automatically by seq2science (v.1.2.4)65 using the atac-seq workflow. Paired-end reads were trimmed with fastp (v.0.23.2)66 with the default options. The genome assembly mm10 was downloaded with genomepy (v.0.16.3)67. Reads were aligned with bwa-mem2 (v.2.2.1)68 with the options ‘-M’. Next, duplicate reads were marked with Picard MarkDuplicates (v.3.0.0) (https://broadinstitute.github.io/picard). General alignment statistics were collected by samtools stats (v.1.16)69. Before peak calling, paired-end info from reads was removed with seq2science so that both mates in a pair get used. Peaks were called with macs2 (v.2.2.7)70 with the options ‘–shift -100 –extsize 200 –nomodel –buffer-size 10000’ in BAM mode. The effective genome size was estimated by khmer (v.3.0)71 by taking the number of unique k-mers in the assembly of the same length as the average read length for each sample. The UCSC genome browser72 was used to visualize and inspect alignment. Quality-control metrics were aggregated by MultiQC (v.1.14)73.
For differential accessibility analysis per-sample MACS2 peak sets were provided to DiffBind (v.3.18.0), which constructed a consensus peak set by merging overlapping peaks across samples. Read counts were generated over the consensus peak regions and normalized using DiffBind default settings. Differential chromatin accessibility was assessed using DiffBind with DESeq2 as the statistical framework, and differentially accessible regions were identified at an adjusted P value (FDR) threshold of 0.05. Significant regions were subsequently annotated to nearby genes using ChIPseeker (annotatePeak) (v.1.44.0) with the TxDb.Mmusculus.UCSC.mm10.knownGene (v.3.10.0) transcript database, using a ±3 kb window around the transcription start sites. Gene identifiers were mapped to gene symbols using org.Mm.eg.db (v.3.21.0). Annotated regions were classified as promoter-associated peaks or distal/enhancer-like peaks (distal, intergenic or intron), and unique gene symbols from each category were extracted for downstream analyses.
To investigate training-associated regulatory programs at individual loci, TPM-normalized RNA expression data were integrated with chromatin accessibility profiles derived from bulk ATAC–seq. Representative genes from the IL-4–STAT6–PPARγ, efferocytosis and pro-resolving lipid mediator programs were examined. Chromatin accessibility profiles were visualized together with published macrophage PPARγ ChIP–seq datasets obtained through ReMap/UCSC (GEO: GSE63696, GSE92606, GSE111854 and GSE107456) and ENCODE candidate cis-regulatory elements. Genomic visualization and manual inspection were performed using the UCSC Genome Browser.
Integration of LEGs with differential expression and chromatin accessibility. LEGs were extracted from significantly enriched upregulated Hallmark pathways identified by GSEA. Gene symbols from all leading-edge subsets were pooled and duplicate entries were removed to generate a unified LEG set. To compare transcriptional and chromatin-level regulation, overlap analyses were performed between LEGs, DEGs and genes associated with DARs. Overlap analyses were performed in R using gene symbols as the common identifier.
scRNA-seq sample preparation and libraries
Mice were euthanized under CO2 and lung tissues were obtained from the control (naive, n = 3), trained (N. brasiliensis infected day 35, n = 3), control + PR8 (3 d.p.i., n = 3) and trained + PR8 (N. brasiliensis trained + PR8 infected, day 38, n = 3) groups. Lungs were injected with 1 ml of collagenase solution (as mentioned above) and incubated for 30 min at 37 °C in an incubator under continuously shaking at 170 rpm. The lungs were crushed using a 70 μm filter and lysed with red blood lysis buffer, then washed with 1× PBS substituted with 2% FBS and 2 mM EDTA. Viable cells were counted for staining with Trypan Blue exclusion using the Countess II (Life Technologies, AMQAX1000) system and stained for viability and lineage markers (CD3, B220, NK1.1), all conjugated to PE-Dazzle, and CD45-BV510. Equal numbers of stained cells from each animal were tagged with lipid oligo tags from the Chromium Next GEM Single Cell 3′ v3.1: Cell Multiplexing kit (10x, 1000261). Cells were pooled into a single sample and sorted using Moflo XDP for scRNA-seq. Cells were sorted using MoFlo cytometer (Beckman Coulter, VXX), and sorted cells were collected in 10% FBS supplemented sterile 1× PBS. A total of 10,500 cells were loaded per channel with an expected yield of 6,000 cells. scRNA-seq libraries were prepared using the Single Cell 3’ Reagent Kits v2 (Chromium Single Cell 3’ Library and Gel Bead Kit v2, PN-120237), Single Cell 3’ Chip Kit v2 PN-120236 and i7 Multiplex Kit PN-120262″ (10x Genomics)3, and the Single Cell 3’ Reagent Kits v2 user guide (manual part no. CG00052 Rev A). Libraries were run on the NovaSeq 6000 system. For alignment, barcode assignment and UMI counting, the Cell Ranger Single Cell Software Suite (v.1.3) was used to perform sample demultiplexing, barcode and unique molecular identifier74 processing and single-cell 3′ gene counting.
Mouse scRNA-seq and normalization
The single-cell count matrix was exported from 10x Genomics cell ranger v.3.0.0 output and subsequently analysed in R v.3.4.1 Single Candle. Clustering and differential expression analysis were performed using the standard functions from Seurat v.2.3.4. Expression markers for each cluster were considered significant with an FDR of 5% and log2[FC] as specified in Fig. 5.
Mouse scRNA-seq data analysis. After confirming the integrity of the cDNA, quality of the libraries, number of cells sequenced and mean number of reads per cell, as a quality control, we used the CellRanger package to map the reads and generate gene–cell matrices. A quality control was then performed on the cells to calculate the number of genes, UMIs and the proportion of mitochondrial genes for each cell; the cells with a low number of covered genes (gene count < 500) and high mitochondrial counts (mt-genes > 0.2) were filtered out. The matrix was then normalized based on library size. A general statistical test was then performed to calculate gene dispersion, base mean and cell coverage to build a gene model for performing PCA. Genes with high coverage (top 500) and high dispersion (dispersion > 1.5) were chosen (1,990 genes) for PCA using iCellR R package (v.1.5.5) (https://CRAN.R-project.org/package=iCellR). t-Distributed stochastic neighbour embedding, UMAP and k-nearest-neighbour-based Network graph drawing Layout (KNetL) were performed. KNetL map has a zoom option that allows the users to see more or less details (more or less subpopulations in cell communities). The network layout used in KNetL map is a force-based layout75 and the zoom option is for changing the force in the system. Force-directed graph drawing algorithms assign attractive (analogous to spring force) and repulsive forces (usually described as analogous to the forces in atomic particles) to separate all pairs of nodes75. The network analysis that has been used in KNetL map have been long used for single-cell analysis and clustering76; here, the nodes of the network layout are extracted and UMAP was performed to create the final plot, a KNetL map. PhenoGraph76 clustering was then performed on the KNetL map results. The marker genes were then found for each cluster and visualized on heat maps, bar plots and box plots. The marker genes were then used to determine the cell types. The proportion (percentage) of the cell communities in each condition were calculated and pseudotime abstract KNetL maps (PAK map) were generated using iCellR.
Pathway and functional analysis. GSEA, pathway and functional analysis were performed using the GSEA tool63. Before performing pathway-level analysis, all genes tested for differential expression were ranked by the DESeq2 Wald statistic. The MSigDB Hallmark gene sets for Mus musculus were then obtained using the msigdbr package (v.25.1.1). Finally, the clusterProfiler package (v.4.16.0) was used to perform the GSEA analysis, whereby the preranked list of genes tested for differential expression was analysed against the MSigDB mice hallmark pathways. Gene sets of size 10–500 were included, and the results were filtered using an FDR threshold of 0.25.
Data availability for raw and analysed data. All the raw ATAC-seq data are deposited into Gene Expression Omnibus (GEO) with the accession number (GSE280183), and all bulk RNA-seq and scRNA-seq data can be found under GEO accession no. GSE343897.
Human scRNA-seq analysis
Human scRNA-seq data were obtained from the publicly accessible IPF Cell Atlas (GEO: GSE136831), comprising lung tissue biopsies from control donors and patients with COPD and IPF. Gene annotations (Ensembl gene IDs and symbols) were obtained using Ensembl BioMart (Ensembl Release 115, September 2025), querying the Human genes (GRCh38.p14) dataset and retrieving the ‘Gene stable ID’ and ‘Gene name’ attributes.
All analyses were performed in R (v.4.4.2) using Seurat (v.5.4.0), edgeR (v.4.4.2), limma (v.3.62.2), fgsea (v.1.32.4), clusterProfiler (v.4.14.6), msigdbr (v.25.1.1), enrichplot (v.1.26.6), dplyr (v.1.1.4), tidyr (v.1.3.2), tibble (v.3.2.0), ggplot2 (v.3.5.1), ggpubr (v.0.6.2), ggrepel (v.0.9.6), patchwork (v.1.3.2) and RColorBrewer (v.1.1.3). Data were subsetted to include only myeloid-annotated cells. The myeloid subset was processed in R with Seurat’s standard workflow, applying NormalizeData(), FindVariableFeatures (vst; n = 2,000), ScaleData() with regression of the Library_Identity covariate, PCA(), FindNeighbors(), FindClusters(), FindSubCluster() and RunUMAP(). Cluster marker genes were identified using Seurat’s FindAllMarkers (), restricting results to positively enriched markers in at least 25% of cells exhibiting a log2[FC] > 0.25; P values were corrected for multiple testing using the Benjamini–Hochberg method, and significance was defined as Padj < 0.05.
Per-individual NAM frequency was calculated as 100 × (NAMs/total myeloid cells). Differences in NAM abundance among the control, COPD and IPF groups were evaluated using a Kruskal–Wallis test followed by pairwise Wilcoxon rank-sum tests with Benjamini–Hochberg FDR correction. To compare NAM transcription across disease states and avoid pseudoreplication from cell-level testing, we performed individual-aware differential expression analysis. NAM cells were aggregated per individual using Seurat’s AggregateExpression(), yielding a gene-by-subject count matrix aligned to metadata (Subject_Identity, Disease_Identity). Low-expressed genes were filtered before modelling, and library-size effects were accounted for using the limma–voom workflow to obtain log2[CPM] values with precision weights. A linear model was fit with a design matrix encoding control, COPD and IPF groups, and contrasts were specified for IPF versus control, COPD versus control and IPF versus COPD. Empirical Bayes moderation was applied, and multiple-testing correction was performed using the Benjamini–Hochberg method, with differentially expressed genes defined at FDR < 0.05. PCA was performed at the individual-level pseudobulk (voom-transformed log2[CPM]) expression matrix; genes were centred and scaled to unit variance before PCA, and variance explained by PC1–PC2 is reported. Samples were coloured by disease_identity (control, COPD, IPF) and plotted in PC1–PC2 space with axis labels indicating the percentage of variance explained by each component; axis limits were fixed to the full range of the data to ensure comparable visualization across panels. Group centroids were computed as the mean PC1 and PC2 coordinates for each disease category and visualized as diamonds to summarize separation among groups.
GSEA of Gene Ontology Biological Process gene sets was performed on the individual-aware differential expression results for NAMs across the three contrasts (IPF versus control, COPD versus control, IPF versus COPD). For each contrast, genes were preranked by the descending log2[FC], with genes lacking statistics removed and, where multiple entries mapped to the same symbol, the instance with the largest absolute log2[FC] was retained. Enrichment was assessed with Benjamini–Hochberg adjustment, with a P-value cut-off of 0.05, and gene set size bounds of 10–500. Moreover, enrichment of our mouse NAM trained gene set was evaluated across the three contrasts. Cross-species gene symbol mapping was performed using Ensembl BioMart (Ensembl Release 115, September 2025), querying the Mouse genes (GRCm39) dataset and retrieving Human Orthologues attributes corresponding to human GRCh38.p14. After conversion of the mouse trained NAM genes to their human orthologues (human gene name; HGNC symbols), fgsea was run with 10,000 permutations and gene set size bounds of 10–500; significance was assessed using Benjamini–Hochberg-adjusted P values (Padj < 0.05).
Thresholds for patients were determined as such: IPF: median, 1.74% (IQR, 1.23%; range, 0.56–6.16%), n = 32; control: median, 0.62% (IQR, 0.74%; range, 0.00–3.87%), n = 28; (Mann–Whitney U-test, FDR-adjusted P = 4.7 × 10−5). COPD: median, 2.16% (IQR, 1.39%; range, 0.00–7.46%), n = 18; versus control: median, 0.62% (IQR, 0.74%; range, 0.00–3.87%), n = 28; (Mann–Whitney U-test, FDR-adjusted P = 3.1 × 10−3) (Fig. 5k). No significant differences were observed between IPF and COPD (FDR-adjusted P = 0.624).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
