Integrated Visualization of Neutrophil Recruitment and Inflammatory Cytokine Profiling in Bronchoalveolar Lavage Fluid

Summary

This study presents an integrated visualization workflow combining standardized bronchoalveolar lavage fluid collection and processing, a seven-marker flow cytometry panel, a 13-analyte multiplex cytokine assay, and a unified analysis pipeline.

Abstract

Bronchoalveolar lavage fluid (BALF) is the most directly accessible sampling matrix for studying lower airway inflammation, yet conventional analyses rely on single-marker flow cytometry combined with analyte-by-analyte enzyme-linked immunosorbent assay (ELISA), an approach that limits information density and obscures cell–factor coupling. This study presents an integrated visualization workflow that couples standardized collection and preprocessing with a seven-marker neutrophil and monocyte flow cytometry panel, a 13-analyte multiplex cytokine assay, and a unified analysis pipeline based on uniform manifold approximation and projection (UMAP), FlowSOM metaclustering, hierarchical clustering, Spearman correlation networks, and principal component analysis (PCA).

The workflow was applied to a single-center prospective cohort of 112 adults undergoing bronchoscopy, comprising 64 patients with microbiologically confirmed pulmonary infection and 48 non-infection controls. Cytological and flow cytometric readouts revealed a pronounced shift from a macrophage-dominant to a neutrophil-dominant alveolar profile during infection, with the neutrophil percentage rising from a median of 3.2% to 68.4% and the absolute neutrophil concentration increasing approximately 50-fold. Multiplex profiling resolved a coordinated three-module cytokine response, with chemokine, proinflammatory, and regulatory modules all upregulated and CXCL8/IL-8, IL-1β, IL-6, CXCL1/GRO-α, and G-CSF showing the strongest single-marker discrimination within the study cohort (area under the curve [AUC] 0.89 to 0.94). Integrative visualization using PCA, a correlation matrix, and a network graph revealed a tightly coupled neutrophil–chemokine axis at the participant level, with BALF neutrophil count emerging as the central network hub.

The workflow uses only material routinely obtained from clinically indicated BAL, requires no additional sample volume, and provides an information-dense framework for resolving airway inflammation in respiratory infection that may be feasible in tertiary-care laboratories equipped with multiparameter flow cytometry and multiplex bead-array platforms, extending into longitudinal, single-cell, and therapeutic-evaluation studies. Broader clinical adoption of individual analytes as biomarkers will require prospective validation in independent multicenter cohorts.

Introduction

Infectious deaths from lower respiratory tract infections remain a significant cause of global mortality, although advancements in vaccine strategies and treatment modalities have been made; nevertheless, this remains especially true among the pediatric and elderly populations¹. Neutrophils are the primary cells reported in the early stages of an innate immune response in the affected respiratory system, a process that is vital for eliminating the invading pathogen; however, uncontrolled neutrophilic inflammation results in the breakdown of the alveolar capillary barrier, oxidative damage to tissue, and the development of acute respiratory distress syndrome (ARDS)^2,3. Neutrophils in the lung airspace also exhibit recruitment kinetics, retention patterns, and effector behavior that differ from other tissues, and these differences influence whether the infection resolves or progresses to severe pneumonia⁴.

The bronchoalveolar lavage fluid (BALF) is an important source of information on the alveolar microenvironment and offers the closest representation of clinical immunity at this site⁵. The chemokine CXCL8 plays an essential role in recruiting neutrophils to the lungs via its receptors, CXCR1 and CXCR2, alongside other chemokines, including CXCL1, CXCL2, and CXCL5⁶. The cytokines IL-1β, IL-6, and TNF-α further enhance this process, whereas IL-10 acts as a counteracting cytokine. Their levels in BALF have been shown to correlate with the severity of lung infections and the development of lung complications due to bacterial, viral, and atypical pneumonia⁷, and multiplex BALF biomarker profiles are being adopted for earlier diagnosis and pathogen discrimination in lower respiratory tract infection⁸.

Routine BALF testing usually fails to fully encompass the biology underlying the test. Cell counts can be estimated by manually differentiating cells in cytospin samples, whereas soluble mediators are measured using individual enzyme-linked immunosorbent assays (ELISAs) and presented as a bar chart. Chemokines, cytokines, and each cell type are evaluated separately, with no correlation established between a particular chemokine concentration gradient and the cell populations that it affects in the sample. The combination of cytokine modules, chemokine-neutrophil dose-response relationships, and patient subpopulations with their own inflammation signatures is often overlooked in the absence of multivariate analysis⁹.

Advances in computation make such a holistic approach progressively realistic. Combining multiparameter flow cytometry and spectral flow cytometry with approaches for dimensionality reduction, such as t-distributed stochastic neighbor embedding (t-SNE) or uniform manifold approximation and projection (UMAP), allows the mapping of high-dimensional immune phenotypes into two-dimensional space, with UMAP proving to have better retention of global topology and reproducibility when applied to cytometry data¹⁰. Multi-analyte panels of cytokines, followed by hierarchical clustering and correlation networks, help identify severity-correlated signatures not discernible with single-cytokine approaches¹¹. Despite progress in the field, no workflow yet brings these methods together to study neutrophil recruitment and cytokine signatures from the same BALF specimens in respiratory infection. This study describes an integrated visualization workflow to investigate neutrophil recruitment using cytology and multiparameter flow cytometry, followed by UMAP and cytokine signature analyses via multiplex immunoassay, with hierarchical clustering, correlation matrices, and principal component analysis. The workflow yields an integrated visualization that encompasses cellular and soluble immune readouts of the airways in a single graph, using BALF samples routinely obtained during clinically indicated bronchoscopy in lower respiratory tract infection.

Protocol

Study population and ethics

This was a single-center prospective observational study conducted at a tertiary teaching hospital between January 2026 and April 2026. The protocol was approved under approval number 2026HL-038, and written informed consent was obtained from every participant before any study procedure.

Adults aged 18 to 75 years undergoing diagnostic flexible bronchoscopy with bronchoalveolar lavage (BAL) as part of routine clinical care were recruited consecutively and assigned to an infection group or a non-infection control group based on their final clinical and microbiological diagnoses. An overview of the integrated visualization workflow, from patient recruitment through to the integrated visual outputs, is shown in Figure 1.

Figure 1: Overview of the integrated visualization workflow for BALF-based profiling of neutrophil recruitment and inflammatory cytokine responses. The five stages summarize, from left to right, prospective recruitment of the infection and control groups with bronchoscopy and BAL, BALF preprocessing into a cell pellet and a supernatant fraction, parallel cellular (Wright–Giemsa cytology and multiparameter flow cytometry) and soluble (13-plex bead-based immunoassay) measurements, the unified visualization pipeline (UMAP with FlowSOM meta-clustering, hierarchical heatmap, principal component analysis, and Spearman correlation network), and the integrated visual output. Abbreviations: BAL = bronchoalveolar lavage; BALF = BAL fluid; UMAP = uniform manifold approximation and projection. Please click here to view a larger version of this figure.

Infection group inclusion criteria were as follows: (1) patients diagnosed with either community-acquired pneumonia (CAP), hospital-acquired pneumonia (HAP) or ventilator-associated pneumonia (VAP) as per Infectious Diseases Society of America (IDSA) and American Thoracic Society (ATS) guidelines; (2) detection of one or more pathogens through BALF culture ≥10⁴ colony-forming units (CFU)/mL, respiratory multiplex polymerase chain reaction (PCR), metagenomic next-generation sequencing of BALF, or specific antigen/antibody testing from paired respiratory/serum specimens.

Control group criteria were as follows: (1) noninfectious reasons for bronchoscopy, consisting of the assessment of an indeterminate lung nodule, abnormal radiology results not suggestive of infection, or examination for chronic cough; (2) no acute pulmonary infection demonstrated by negative results from microbiological testing and lack of any signs of inflammation or systemic infection at the time of bronchoscopy; and (3) no use of antibiotics in the 30 days prior to participation. Controls were recruited from patients undergoing bronchoscopy for other clinically indicated pulmonary conditions rather than from healthy volunteers, reflecting the population from which clinically indicated BAL is routinely obtained at the study center.

Exclusion criteria included: (1) malignancy and/or cytotoxic chemotherapy in the preceding six months; (2) HIV seropositivity; (3) solid organ or hematopoietic stem cell transplantation; (4) sustained immunosuppressive therapy involving a daily equivalent prednisone dosage > 0.3 mg/kg for more than 3 weeks; (5) pregnant/lactating females; (6) history of bleeding diathesis or any contraindications to flexible bronchoscopy; (7) pre-existing severe chronic respiratory disease (chronic obstructive pulmonary disease [COPD] Global Initiative for Chronic Obstructive Lung Disease [GOLD] stage III/IV, idiopathic pulmonary fibrosis, severe persistent asthma, or clinically significant bronchiectasis); (8) active tuberculosis; (9) involvement in another interventional study within the last 30 days. BAL was performed during the first bronchoscopy, within 72 h of admission for infections, and within 48 h of intubation in mechanically ventilated patients with infections.

Sample size calculations were focused on the primary cellular outcome: the number of neutrophils in BALF among groups. With respect to the anticipated effect size of 0.6 SD between infected and non-infected airways, and considering a 15% dropout allowance, a minimum sample size of 45 subjects per group was needed for two-sided α = 0.05 with 80% statistical power. Out of 156 eligible candidates, 27 did not satisfy the eligibility requirements (18 failed to fulfill the inclusion criteria mainly due to inability to identify a pathogenic agent responsible for their disease within the stipulated period or lack of agreement with the procedure itself; 9 breached exclusion criteria primarily owing to active malignancy or history of chemotherapeutic administration, chronic corticosteroid use, severe pre-existing chronic lung disease, or bleeding diathesis). The remaining 129 participants underwent bronchoscopy. Subsequently, 17 subjects were withdrawn from the study (7 because of inadequate BALF yield < 40%; 4 because of consent withdrawal while hospitalized; and 6 because of reclassification of subjects after final microbiologic evaluation – 4 control candidates diagnosed with occult infection based on culture results; and 2 infection candidates who lacked evidence of infection after culture analysis). The final pool of study subjects comprised 112 participants: 64 with infection and 48 without.

BALF collection and preprocessing

Bronchoscopy was performed by skilled operators using the standard institutional protocol, consistent with published technical recommendations for bronchoalveolar lavage¹². Lidocaine was used topically, and patients received moderate intravenous sedation while supplemental oxygen was provided to keep oxygen saturation (SpO₂) >92%. The flexible bronchoscope was inserted and positioned into the bronchopulmonary segment with radiological abnormalities. In cases where the radiological abnormalities were diffuse or unlocalizable, right middle lobe or left lingular segments were chosen. After each instillation, gentle manual suction with negative pressure < 100 mmHg was performed. A total of three to five instillations of sterile saline at 37 °C, 20–40 mL each, for a total volume of 100–150 mL, were administered. The recovery ratio was estimated as the recovered volume relative to the instilled volume, and only samples with a recovery ≥40% were considered for analysis.

The collected samples from the second and subsequent washes were pooled to reduce bronchial contamination and obtain an alveolar-enriched sample. The bronchoalveolar lavage fluid (BALF) pool was kept on ice and processed within 2 h of collection, taking into account recent findings on changes in the cellular composition and differential counts after prolonged storage and warming of human BALF samples¹³. After filtration using a sterile nylon mesh with a pore size of 70 microns to filter out the mucus, performed by mechanical passage alone without an added mucolytic or anticlumping agent such as dithiothreitol, since pilot testing showed this was sufficient to remove mucus clumps without measurable loss of downstream surface-epitope integrity, the fluid was subjected to centrifugation at 400 × g for 10 min at 4 °C. The supernatant was decanted and aliquoted into low-protein-binding cryovials for cytokine analysis, and the cells were suspended in phosphate-buffered saline (PBS) with 0.5% bovine serum albumin and 2 mM ethylenediaminetetraacetic acid (EDTA), with part of them kept for cytology and flow cytometry. Cytokine samples were frozen in liquid nitrogen, kept at –80 °C, and subjected to no more than one freeze-thaw cycle before use.

Neutrophil cytology and flow cytometric quantification

In terms of cytopathology, 100 µL of the BALF cell suspension was applied to the glass slide using a cytocentrifuge (see Table of Materials) at 300 × g for 5 min, air-dried, and stained with Wright-Giemsa. Differential counts were performed by a hematology technologist without knowledge of the patient's medical background, who counted a minimum of 400 cells per slide and calculated percentages of macrophages, lymphocytes, neutrophils, and eosinophils. BALF neutrophil percentages were considered the primary outcome measure for cytopathology, given current evidence suggesting their usefulness in distinguishing infectious from noninfectious lung diseases in patients with pulmonary infiltrates¹⁴.

The flow cytometry quantification protocol was based on a recent validation of the four-color bronchoalveolar lavage fluid (BALF) leukocyte differentiation protocol¹⁵ and included additional markers to identify neutrophils. Specifically, around 5 × 10⁵ cells isolated from BALF were blocked with a human Fc receptor-blocking reagent (see Table of Materials) for 10 min at 4 °C, and then stained for 20 min at 4 °C with antibodies against CD45 (clone HI30), CD66b (clone G10F5), CD15 (clone HI98), CD16 (clone 3G8), CD11b (clone ICRF44), CD14 (clone M5E2), and HLA-DR (clone L243), in addition to a fixable viability dye (see Table of Materials). For erythrocyte lysis, a commercial erythrocyte lysing solution (see Table of Materials) was used for 10 min at room temperature. After erythrocyte lysis, cells were washed twice with PBS supplemented with 2% fetal calf serum and resuspended in flow buffer, then analyzed with a multicolor flow cytometer (see Table of Materials). Compensation was performed using single-stained compensation beads (see Table of Materials), and fluorescence-minus-one and unstained tubes were included for each panel.

At least 1 × 10⁵ CD45⁺ events were acquired for each sample. Singlets on the FSC-A/FSC-H plot were gated first, followed by live (viability dye⁻) CD45⁺ leukocytes, after which exclusion of CD14⁺HLA-DR⁺ monocyte-macrophages was performed to select high-side scatter CD66b⁺ granulocytes. Neutrophils were then defined within the granulocyte population as CD15⁺CD16⁺CD11b⁺ events. This primary CD66b⁺ granulocyte gate was drawn broadly around all high-side-scatter events rather than being restricted to the CD15⁺CD16⁺CD11b⁺ boundary, so that CD16-dim, CD11b-high activated neutrophils and other less mature granulocyte forms remained within the total automated neutrophil count; the stricter CD15⁺CD16⁺CD11b⁺ definition was applied only to characterize the mature-neutrophil subset reported separately in the downstream clustering analyses. The neutrophil count in BALF (cells/µL) was estimated based on neutrophil percentage and the total nucleated cell count determined using an automated hematology analyzer (see Table of Materials). Compensation matrices, gating strategy, and individual sample quality-control reports were archived, and analysts remained blind to participant group assignments throughout acquisition and analysis.

Multiplex inflammatory cytokine profiling

BALF supernatants were thawed once on wet ice, followed by brief vortexing and clarification by centrifugation at 10,000 × g for 5 min at 4 °C directly prior to assay. Measurement of 13 analytes in three functional categories was achieved using a multiplex bead-based immunoassay platform with a custom magnetic bead panel (see Table of Materials). The pro-inflammatory group comprised IL-1β, IL-6, IL-17A and TNF-α; the chemokine group comprised CXCL8/IL-8, CXCL1/GRO-α, CXCL10/IP-10, CCL2/MCP-1 and CCL3/MIP-1α; and the regulatory group comprised IL-10, IFN-γ, IL-1RA and G-CSF. The choice of analytes was based on those most commonly reported as informative in multiplex analyses of BALF or serum cytokines in pneumonia/severe respiratory infection studies^7,11. The dilution factor (1:2 in assay buffer) was determined from a pilot experiment conducted on a representative subset of both infected and control samples to ensure that >90% of measurements fell within the standard curve.

The standard curves were determined using seven-point serial dilutions of the manufacturer's recombinant standards, prepared in duplicate on every plate, including one blank well to subtract the background. The curve fitting was conducted via five-parameter logistic (5PL) regression by the instrument's proprietary acquisition and analysis software (see Table of Materials). The lower limit of quantification (LLOQ) was calculated as the smallest amount of standard yielding a back-fit recovery between 80% and 120% and an intrarun coefficient of variation (CV) ≤ 20%. The values that fell below the LLOQ were substituted by LLOQ/√2. Two pooled bronchoalveolar lavage fluid (BALF) quality-control samples, each from a representative infected patient and an unaffected subject, were aliquoted and added to each plate to track the long-term performance. Quality assurance standards were based on best-practice recommendations derived from a 12 year multisite multiplex bead-based assay proficiency program¹⁶, including requirements for intra-assay CV ≤ 10% in duplicate wells, inter-plate CV ≤ 20% for quality control, and a minimum bead count of ≥ 50 per analyte per well. Plates that failed to meet the above criteria were rerun. Analysis was performed in a single batch within four weeks, using a single kit lot, with random interleaving of the infection and control samples.

Visualization analysis pipeline

The Flow Cytometry Standard (FCS) files were first analyzed using commercial flow cytometry software (see Table of Materials) to compensate and gate live CD45⁺ singlets, and then exported for unbiased analysis in R v4.4. For each sample, 10,000 CD45⁺ events were downsampled and concatenated, and signal intensities were transformed using the arcsinh function with a cofactor of 150. Dimension reduction and UMAP projections were conducted using the uwot package at 15 nearest neighbors and a minimum distance of 0.1, along with parallel t-SNE embeddings using the Rtsne package with perplexity 30¹⁰. Phenotype clusters were identified using FlowSOM with 100 self-organizing map nodes and further divided into 12–15 clusters via meta-clustering, using an optimization-based approach that has recently been validated for high-dimensional cytometry data¹⁷. Each cluster was annotated based on its median expression of lineage and activation markers and reviewed against the manual gates.

For the analysis of multiplex cytokine data, raw values were first log₁₀ transformed and then standardized between analytes via z-scoring before additional analyses. Both analytes and patients underwent hierarchical clustering using Euclidean distance and Ward.D2 linkage; the results were displayed in a row- and column-clustered heatmap. Such an approach was used to detect co-expressed cytokine modules and classify patient subsets in multiplex inflammatory cohorts⁹. Correlation matrices were generated for all cellular and soluble biomarkers using pairwise Spearman rank correlations. Heatmaps of the correlation matrices were plotted alongside a network in which nodes correspond to biomarkers and edges indicate correlations that survived multiple-testing correction using the Benjamini–Hochberg (BH) method. Edge width was scaled according to the correlation coefficient magnitude |ρ|, and edges were included only if |ρ| ≥ 0.30 and p < 0.05 after correction. Principal component analysis (PCA) was performed on the standardized matrix containing cells and cytokines using the prcomp() function with unit variance scaling. Participants were projected onto the first two principal components with infection and control groups overlaid. Consistent downsampling and standardization of all analyses ensured comparability between cellular and soluble biomarker readouts across all visualizations.

Statistical analysis

Continuous variables were described using the median (interquartile range, IQR) when non-normal according to Shapiro–Wilk tests, and mean (standard deviation, SD) otherwise. Categorical variables were described using frequencies and proportions. Comparisons between groups were performed using the Mann–Whitney U test when continuous variables had non-normal distributions, the Welch t-test when continuous variables were normally distributed, and the chi-squared test/Fisher's exact test when applicable for categorical variables. The associations between continuous variables were assessed through Spearman's rank correlation coefficient. To account for the false discovery rate (FDR) in multiplex cytokine and correlation analyses, p-values were adjusted using the Benjamini–Hochberg method with a cutoff of q < 0.05. Two-sided p < 0.05 was considered statistically significant in all primary analyses. Where possible, cases with missing data on individual variables were retained and analyzed by complete-case analysis for the affected variable, while multiple imputations using chained equations (m = 20) were performed to conduct sensitivity analyses. Receiver operating characteristic curves and the corresponding area under the curve (AUC) were generated for individual cytokines using the pROC package within the full study cohort; a separate external validation cohort was not used. Pre-specified sensitivity analyses repeated the primary group comparisons after excluding participants who had received antibiotics within 72 h before bronchoscopy, and after stratifying the infection group by mechanical ventilation status and by pneumonia subtype (community-acquired, hospital-acquired, or ventilator-associated). Analyses were performed using the rstatix, FlowSOM, uwot, ComplexHeatmap, igraph, pROC, and mice packages in R v4.4. Figures were generated using commercial scientific graphing software (see Table of Materials).

Results

Study population and sample characteristics

There were 112 subjects in the final analysis cohort, with 64 subjects in the infection group and 48 subjects in the non-infection control group. Baseline demographics were similar between groups, while laboratory parameters related to infection were higher in the infection group, as shown in Table 1. Infections were primarily caused by Gram-negative bacteria, with the main bacterial pathogens being Klebsiella pn...

Discussion

The robustness of this protocol depends on a few decisions made before and during analysis. Setting the lower-bound recovery rate at 40% and combining the second and subsequent aliquots was an intentional decision, as the first aliquot contains a higher enrichment of bronchial airway materials, while the recovery rate systematically varies with age, sex, smoking history, and the target bronchus¹⁸. Blocking Fc receptors and titrating antibodies against individual clones to ensure correct gating bou...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Anti-human CD11b antibody, Brilliant Violet 421, clone ICRF44	BioLegend	Cat. No. 301323	Neutrophil activation marker; violet laser channel URL: https://www.biolegend.com
Anti-human CD14 antibody, APC/Cyanine7, clone M5E2	BioLegend	Cat. No. 301820	Monocyte-macrophage exclusion marker; red laser channel URL: https://www.biolegend.com
Anti-human CD15 antibody, PE, clone HI98	BioLegend	Cat. No. 301906	Neutrophil lineage marker; blue laser channel URL: https://www.biolegend.com
Anti-human CD16 antibody, Brilliant Violet 510, clone 3G8	BioLegend	Cat. No. 302048	Neutrophil maturation and NK cell marker; violet laser channel URL: https://www.biolegend.com
Anti-human CD45 antibody, PerCP, clone HI30	BioLegend	Cat. No. 304026 / RRID: AB_893338	Pan-leukocyte marker; blue laser channel URL: https://antibodyregistry.org/AB_893338
Anti-human CD66b antibody, FITC, clone G10F5	BioLegend	Cat. No. 305103	Granulocyte marker; blue laser channel URL: https://www.biolegend.com
Anti-human HLA-DR antibody, PerCP/Cyanine5.5, clone L243	BioLegend	Cat. No. 307629	Antigen-presenting cell marker; blue laser channel URL: https://www.biolegend.com
Automated haematology analyser (Sysmex XN-1000)	Sysmex Corporation	Instrument (model XN-1000)	Used to determine total nucleated cell count for absolute neutrophil calculation URL: https://www.sysmex.com
Bovine serum albumin	Sigma-Aldrich (MilliporeSigma)	Cat. No. A7906	Used in flow cytometry staining buffer (0.5% w/v in PBS with 2 mM EDTA) URL: https://www.sigmaaldrich.com/US/en/product/sigma/a7906
ComplexHeatmap package (R/Bioconductor)	Bioconductor Project	v2.20 (open-source)	Used for hierarchical clustering heatmaps of multiplex cytokine data URL: https://bioconductor.org/packages/ComplexHeatmap
Cytocentrifuge (Cytospin 4)	Thermo Fisher Scientific	Cat. No. A78300003	Used to prepare cytology slides from BALF cell suspension URL: https://www.thermofisher.com/order/catalog/product/A78300003
Erythrocyte lysing solution (BD FACS Lysing Solution, 10X concentrate)	BD Biosciences	Cat. No. 349202	Used for erythrocyte lysis after surface staining URL: https://www.bdbiosciences.com
Ethylenediaminetetraacetic acid (EDTA), 0.5 M pH 8.0	Thermo Fisher Scientific (Invitrogen)	Cat. No. AM9260G	Used in flow cytometry staining buffer at 2 mM URL: https://www.thermofisher.com/order/catalog/product/AM9260G
Fixable viability dye (Zombie NIR Fixable Viability Kit)	BioLegend	Cat. No. 423105	Live-dead exclusion prior to gating; red laser channel URL: https://www.biolegend.com
Flexible video bronchoscope	Olympus Corporation (BF-1TQ290)	Instrument (model BF-1TQ290)	Used for BAL under moderate sedation URL: https://www.olympus-global.com
Flow cytometer, multicolour (CytoFLEX LX)	Beckman Coulter Life Sciences	Instrument (model CytoFLEX LX, six-laser configuration)	Used for acquisition of the seven-marker neutrophil/monocyte flow panel URL: https://www.beckman.com
Flow cytometry analysis software (FlowJo)	BD Biosciences (FlowJo LLC)	v10.10	Used for compensation and manual gating of FCS files URL: https://www.flowjo.com
FlowSOM package (R/Bioconductor)	Bioconductor Project	v2.12 (open-source)	Used for self-organising-map clustering and meta-clustering of flow cytometry data URL: https://bioconductor.org/packages/FlowSOM
GraphPad Prism (statistical graphing software)	GraphPad Software (Dotmatics)	v10	Used to generate final figures URL: https://www.graphpad.com
Human Fc receptor-blocking reagent (Human TruStain FcX)	BioLegend	Cat. No. 422301 / RRID: AB_2818986	Fc receptor block prior to antibody staining URL: https://antibodyregistry.org/AB_2818986
igraph package (R)	CRAN	v2.0 (open-source)	Used to construct and visualise the correlation network graph URL: https://cran.r-project.org/package=igraph
Low-protein-binding cryogenic vials, 2.0 mL	Thermo Fisher Scientific (Nunc)	Cat. No. 375418	Used for BALF supernatant aliquots for cytokine analysis URL: https://www.thermofisher.com/order/catalog/product/375418
Luminex multiplex bead-based immunoassay platform (MAGPIX)	Luminex Corporation (DiaSorin)	Instrument (model MAGPIX)	Used for acquisition of the multiplex cytokine assay URL: https://www.luminexcorp.com
May-Grünwald-Giemsa stain solution	Sigma-Aldrich (MilliporeSigma)	Cat. No. GS500	Used for cytospin differential cell count URL: https://www.sigmaaldrich.com/US/en/product/sigma/gs500
mice package (R)	CRAN	v3.16 (open-source)	Used for multiple imputation by chained equations in sensitivity analysis URL: https://cran.r-project.org/package=mice
Multiplex magnetic bead cytokine/chemokine panel (13-plex, custom; MILLIPLEX MAP Human Cytokine/Chemokine)	Merck / MilliporeSigma	Cat. No. HCYTOMAG-60K	Custom panel measuring the 13 analytes described in Section 2.4 URL: https://www.merckmillipore.com
Nylon mesh cell strainer, 70 µm	Corning (Falcon)	Cat. No. 352350	Used for mechanical filtration of BALF to remove mucus clumps URL: https://www.corning.com
Phosphate-buffered saline (PBS), pH 7.4	Thermo Fisher Scientific (Gibco)	Cat. No. 10010023	Base buffer for cell resuspension and staining URL: https://www.thermofisher.com/order/catalog/product/10010023
pROC package (R)	CRAN	v1.18 (open-source)	Used to generate ROC curves and AUC values URL: https://cran.r-project.org/package=pROC
R (statistical computing environment)	R Core Team	v4.4 (open-source)	Platform for statistical analyses and visualisation URL: https://www.r-project.org
rstatix package (R)	CRAN	v0.7 (open-source)	Used for group comparisons and correlation statistics URL: https://cran.r-project.org/package=rstatix
Rtsne package (R)	CRAN	v0.17 (open-source)	Used for t-SNE embeddings URL: https://cran.r-project.org/package=Rtsne
Single-stained compensation beads (UltraComp eBeads Compensation Beads)	Thermo Fisher Scientific (Invitrogen)	Cat. No. 01-2222-42	Used for single-colour compensation controls URL: https://www.thermofisher.com/order/catalog/product/01-2222-42
Sodium chloride solution for injection, 0.9% (sterile saline)	Baxter International	Cat. No. 2B1324	Used for BAL instillation (100–150 mL, 20–40 mL aliquots at 37 °C) URL: https://www.baxter.com
uwot package (R)	CRAN	v0.2 (open-source)	Used for UMAP dimensionality reduction URL: https://cran.r-project.org/package=uwot
xPONENT (multiplex assay acquisition and curve fitting software)	Luminex Corporation (DiaSorin)	v4.3	Instrument-proprietary software for 5PL standard curve fitting URL: https://www.luminexcorp.com