Pathobiological signatures of dysbiotic lung injury in pediatric patients undergoing stem cell transplantation

Abstract

Hematopoietic cell transplantation (HCT) uses cytotoxic chemotherapy and/or radiation followed by intravenous infusion of stem cells to cure malignancies, bone marrow failure and inborn errors of immunity, hemoglobin and metabolism. Lung injury is a known complication of the process, due in part to disruption in the pulmonary microenvironment by insults such as infection, alloreactive inflammation and cellular toxicity. How microorganisms, immunity and the respiratory epithelium interact to contribute to lung injury is uncertain, limiting the development of prevention and treatment strategies. Here we used 278 bronchoalveolar lavage (BAL) fluid samples to study the lung microenvironment in 229 pediatric patients who have undergone HCT treated at 32 children's hospitals between 2014 and 2022. By leveraging paired microbiome and human gene expression data, we identified high-risk BAL compositions associated with in-hospital mortality (P = 0.007). Disadvantageous profiles included bacterial overgrowth with neutrophilic inflammation, microbiome contraction with epithelial fibroproliferation and profound commensal depletion with viral and staphylococcal enrichment, lymphocytic activation and cellular injury, and were replicated in an independent cohort from the Netherlands (P = 0.022). In addition, a broad array of previously occult pathogens was identified, as well as a strong link between antibiotic exposure, commensal bacterial depletion and enrichment of viruses and fungi. Together these lung-immune system-microorganism interactions clarify the important drivers of fatal lung injury in pediatric patients who have undergone HCT. Further investigation is needed to determine how personalized interpretation of heterogeneous pulmonary microenvironments may be used to improve pediatric HCT outcomes.

Fig. 1. Study design and clinical outcomes.

a, Patients were recruited from 32 participating children’s hospitals in the United States, Canada and Australia. b, Study design diagram. c, BAL processing and analysis workflow. d, Four microbiome–transcriptome clusters were identified. e, In-hospital survival for all patients (left) and the subset requiring respiratory support before testing (right) was plotted according to BAL cluster; differences were analyzed with the log-rank test.

Fig. 2. The BAL microbiome.

a, The fraction (left) and mass (right) of major bacterial, viral and fungal phyla were plotted, with the shading representing the average for each of the four BAL clusters (n = 127, 74, 45 and 32 for clusters 1–4, respectively). The average mass of bacterial genera and species in each of the four BAL clusters are shown on the right. b, Taxonomic richness and diversity were plotted across the four BAL clusters. Richness and diversity varied across clusters (Kruskal–Wallis test, P < 0.001 and P = 0.002, respectively). c, Microorganisms associated with in-hospital mortality were identified using negative binomial generalized linear models (edgeR R package) and were plotted according to the log fold change (position, color) and FDR (dot size). d, Taxonomic richness and Simpson’s alpha diversity stratified according to survival status at the time of the most recent BAL (n = 184 survivors, n = 45 nonsurvivors). Richness and diversity differed according to survival outcome (Wilcoxon rank-sum test, P = 0.025 and P = 0.006, respectively). e, Microbial alignments to the KEGG metabolic pathways were averaged for each BAL cluster. f, Selected metabolic pathways that differed across the BAL clusters are shown. log₁₀-normalized expression varied across clusters (Kruskal–Wallis test, FDR < 0.001 for each of glycolysis/gluconeogenesis, oxidative phosphorylation, fatty acid biosynthesis and butanoate metabolism). For all box plots: the boxes indicate the median and IQR; the whiskers extend to the largest value above the 75th percentile (or smallest value below the 25th percentile), that is, within 1.5 times the IQR.

Fig. 3. BAL pathogen detection.

a, Left: dot plots of common community-transmitted respiratory viruses (left), herpesviruses (middle) and all other viruses (right) detected in the cohort, plotted according to microbial mass (x axis) and microbiome dominance (y axis). Right: bar chart comparing viral detection across the four BAL clusters according to hospital tests and metagenomic sequencing. b, Left: all H. influenzae, S. aureus and S. pneumoniae detected in the cohort were plotted, with the dashed lines indicating the cutoffs of mass ≥10 pg and bacterial dominance ≥20%. Taxa above these cutoffs are shown in the upper-right quadrant (shaded in yellow) to indicate outliers within the cohort. Right: bar chart comparing potentially pathogenic bacteria detected across the four BAL clusters according to hospital tests and metagenomic sequencing. c, Left: all microorganisms detected in the BAL of three patients are shown, with the arrows indicating fungi present in high quantities. Right: bar chart comparing potentially pathogenic eukaryotes detected across the four BAL clusters according to hospital tests and metagenomic sequencing.

Fig. 4. Antibiotic exposure and impact on BAL microbiome.

a, Days of antimicrobial exposure are listed for antibacterials (black), antifungals (green) and antivirals (blue). Patients are listed in the columns and the shading indicates the number of days of exposure to each antibiotic in the week preceding BAL. b, AES was calculated before each BAL as the sum of antibiotic exposure days × a broadness weighting factor, summed for all therapies received in the week preceding BAL. AES varied across the clusters (n = 127, 74, 45 and 32 for clusters 1–4, respectively) and was highest for patients in cluster 4 (Kruskal–Wallis test, P = 0.005). For all box plots: the boxes indicate the median and IQR; the whiskers extend to the largest value above the 75th percentile (or the smallest value below the 25th percentile), that is, within 1.5 times the IQR. c, Negative binomial generalized linear models were used to test for BAL microorganisms associated with AES. Microorganisms are listed in the rows, with phyla shown on the left and bacterial genera shown on the right.

Fig. 5. BAL gene expression.

a, DEGs were identified using a four-way analysis of variance-like analysis with negative binomial generalized linear models. Mean normalized expression levels for significant genes are displayed for the four BAL clusters. b, Individual DEGs were identified across the four clusters (edgeR R package); variance-stabilized transformed gene counts for select genes highest in each of the four clusters were plotted (n = 127, 74, 45 and 32 for clusters 1–4, respectively). For all box plots: boxes indicate the median and IQR; the whiskers extend to the largest value above the 75th percentile (or smallest value below the 25th percentile), that is, within 1.5 times the IQR. c, Gene set enrichment scores to REACTOME pathways were calculated and example gene sets most enriched in each of the four clusters are shown.

Extended Data Fig. 1. Microbial KEGG Metabolism Pathways.

Mean ERCC-transformed normalized KEGG pathway expression for microbial Carbohydrate, Energy, Lipid, and Glycan metabolism pathways.

Extended Data Fig. 2. Antimicrobial Resistance Gene Expression By BAL Cluster.

Antimicrobial resistance gene (AMR) expression was derived from human-subtracted sequencing files processed using the. CZID Antimicrobial Resistance (AMR) Gene Pipeline v0.2.4-beta, which leverages the Resistance Gene Identifier (RGI) v6.0.0 to generate read k-mer alignments (KMA) against the Comprehensive Antibiotic Resistance Database (CARD) v3.2.3 and WILDCARD v3.1.0. AMR transcripts were summed across all AMR genes and normalized to sample ERCC reads (left) and additionally to total BAL mass of bacteria (right) and varied by cluster (Kruskal-Wallis p < 0.001 and p < 0.001, respectively). n = 127, 74, 45, and 32 for Clusters 1-4, respectively. For all box-whisker plots: boxes indicate the median and interquartile range and whiskers extend to the largest value above the 75th percentile (or smallest value below the 25th percentile) that is within 1.5 times the IQR.

Extended Data Fig. 3. In-Hospital Survival Stratified by BAL Pathogen Detected.

In-hospital survival for patients with a viral pathogen (top), bacterial pathogen (middle), or eukaryotic pathogen (bottom) detected on BAL, relative to no pathogen detected on BAL. Survival curves were compared using the two-sided log-rank test.

Extended Data Fig. 4. In-Hospital Survival Stratified by Antibacterial Exposure Score.

Antibacterial exposure score (AES) was divided into 4 quartiles of equal patient number and in-hospital survival was plotted for each quartile and compared with the two-sided log-rank test.

Extended Data Fig. 5. BAL Cell-Type Deconvolution and Imputed Cell-Specific Gene Expression.

Cell type fractions were imputed using bulk gene expression and a reference single cell lung atlas. (a) Mean centered, scaled fractions are depicted for each cluster and varied across clusters for Monocyte/macrophages, neutrophils, CD4 + T-cells, and CD8+ Tcells (Kruskal-Wallis p = 4.14 × 10-24, 3.81×10-5, 1.89×10-22, and 3.36×10-13, respectively). (b) Raw values are shown for specific cell types. (c) Monocyte/macrophage-specific expression of the ‘GOBP Myeloid Leukocyte Activation’ gene set was imputed. Genes that were statistically significantly differentially expressed across clusters were selected for the heatmap, and average cell-type specific gene expression across the 4 clusters is displayed. (d) Lymphocyte-specific expression of ‘GOBP Lymphocyte Activation’ gene set was imputed. Genes that were statistically significantly differentially expressed across clusters were selected for the heatmap, and average cell-type specific gene expression across the 4 clusters is displayed. n = 127, 74, 45, and 32 for Clusters 1-4, respectively. For all box-whisker plots: boxes indicate the median and interquartile range and whiskers extend to the largest value above the 75th percentile (or smallest value below the 25th percentile) that is within 1.5 times the IQR.

Extended Data Fig. 6. BAL T-Cell Receptor Repertoires.

CDR3 alignments were computed and clonotypes and Shannon diversity of TCRα alignments are shown for each of the BAL clusters. n = 127, 74, 45, and 32 for Clusters 1-4, respectively. TRA clonotypes and shannon diversity varied by cluster (Kruskal-Wallis p = 1.32 × 10-7 and 3.27 × 10-8, respectively). For all box-whisker plots: boxes indicate the median and interquartile range and whiskers extend to the largest value above the 75th percentile (or smallest value below the 25th percentile) that is within 1.5 times the IQR.

Extended Data Fig. 7. BAL Cluster Transitions.

34 patients had ≥2 BALs in this study (total 49 BALs were repeat samples). Cluster transitions are shown here, indicating a general transition away from the low-risk Cluster 1 on repeat samples.

Extended Data Fig. 8. Validation Cohort Survival Stratified by BAL Cluster.

A random forest classifier using BAL metagenomic and transcriptomic data was grown using the derivation validation set. The classifier was applied to BAL data from a validation cohort, and 1-year non-relapse mortality was plotted according to cluster assignment and compared using the two-sided log-rank test.

Extended Data Fig. 9. BAL Cluster Schema.

(a) BAL Cluster 1 was most common, had moderate microbial burden, low rates of infection, predominantly alveolar macrophage-related signaling, and the lowest mortality rates. (b) Cluster 2 showed high rates of microbial burden and bacterial infections, higher neutrophil markers, and moderate mortality. (c) Cluster 3 showed microbiome depletion with enrichment of viruses and fungi and fibroproliferative gene expression. (d) Cluster 4 showed significant microbiome depletion with relative sparing of Staphylococci and enrichment of viruses, commensurate with lymphocytic inflammation, cellular injury, and the highest mortality rate.