Pulmonary microbiome and gene expression signatures differentiate lung function in pediatric hematopoietic cell transplant candidates

Abstract

Impaired baseline lung function is associated with mortality after pediatric allogeneic hematopoietic cell transplantation (HCT), yet limited knowledge of the molecular pathways that characterize pretransplant lung function has hindered the development of lung-targeted interventions. In this study, we quantified the association between bronchoalveolar lavage (BAL) metatranscriptomes and paired pulmonary function tests performed a median of 1 to 2 weeks before allogeneic HCT in 104 children in The Netherlands. Abnormal pulmonary function was recorded in more than half the cohort, consisted most commonly of restriction and impaired diffusion, and was associated with both all-cause and lung injury-related mortality after HCT. Depletion of commensal supraglottic taxa, such as Haemophilus, and enrichment of nasal and skin taxa, such as Staphylococcus, in the BAL microbiome were associated with worse measures of lung capacity and gas diffusion. In addition, BAL gene expression signatures of alveolar epithelial activation, epithelial-mesenchymal transition, and down-regulated immunity were associated with impaired lung capacity and diffusion, suggesting a postinjury profibrotic response. Detection of microbial depletion and abnormal epithelial gene expression in BAL enhanced the prognostic utility of pre-HCT pulmonary function tests for the outcome of post-HCT mortality. These findings suggest a potentially actionable connection between microbiome depletion, alveolar injury, and pulmonary fibrosis in the pathogenesis of pre-HCT lung dysfunction.

Figure 1.. Association between Pre-HCT PFTs, Clinical Characteristics, and Post-HCT Outcomes.

(A) Box-and-whisker plot of age distribution are shown for patients with restriction (n=35), obstruction (n=12), or neither (n=57) using conventional PFT categorization. Patients with either restriction or obstruction were older than patients with normal PFTs. Data were analyzed using Dunn’s Test. (B) Box-and-whisker plot of PFT distributions are shown according to patient sex (females, n=43; males, n=61). Data were analyzed using Wilcoxon rank-sum tests. Box-and-whisker plots in (A) and (B) depict the median and interquartile range (box) with 1.5 times the interquartile range (whiskers). (C) Kaplan-Meier post-HCT survival estimates are shown according to pre-HCT FVC_%pred (left) and TLC_%pred (right). Both FVC_%pred and TLC_%pred were associated with post-HCT survival. Data were analyzed by Cox regression.

Figure 2.. Association between BAL Microbiome and FVC_%pred.

(A) The volcano plot shows the association between BAL microbe RNA mass and FVC_%pred. Associations were tested using negative binomial generalized linear models accounting for patient age and sex covariates. Positive log₂ fold change (x-axis) indicates that greater microbial RNA mass is associated with greater FVC_%pred; negative log₂ fold change indicates that greater microbial RNA mass is associated with lower FVC_%pred. (B and C) Box-and-whisker plot of variance stabilizing transformed (vst) RNA masses of Haemophilus (B) and Staphylococcus (C) are shown according to category of FVC_%pred (>100%, n=25; 80–100%, n=40; 60–80%, n=28; <60%, n=11). Worse restriction is associated with lower Haemophilus and greater Staphylococcus RNA mass in BAL fluid.Box-and-whisker plots in (B) and (C) depict the median and interquartile range (box) with 1.5 times the interquartile range (whiskers).

Figure 3.. Association between BAL Hallmark Gene Expression and Pulmonary Function.

Hallmark gene set variance analysis using Poisson gene count distributions was used to calculate gene set enrichment scores for each patient. Spearman rank-based correlations between gene set enrichment scores (rows) and PFTs (columns) were calculated and plotted (cells) with hierarchical clustering of columns and rows according to Euclidean distance. Dendrogram-based clustering shows that PFT measures of lung size and diffusion were positively associated with a cluster of hallmark gene sets corresponding to immune activation, cellular growth and division, and metabolism and were negatively associated with a cluster of hallmark gene sets corresponding to mesenchymal transition, cellular fate and differentiation, and the lung stroma.

Figure 4.. Association between BAL Gene Expression and Pulmonary Function Testing.

Heatmaps of differential gene expression are plotted according to pulmonary restriction (left) and impaired diffusion (right). Negative binomial models were used to identify differentially expressed genes (rows) across the continuum of FVC_%pred and DLCO/Va_%pred (columns). Cell shading indicates mean gene expression for each PFT subgroup after variance stabilizing transformation, centering, and scaling. Dendrogram-based clustering shows that declining FVC_%pred was associated with a decrease in expression of immune-related genes and an increase in expression of airway epithelial genes relating to surfactant, keratinization, and olfactory responses. Declining DLCO/Va_%pred was also associated with a decrease in expression of immune-related genes and an increase in expression of airway epithelial genes as well as hormone response genes.

Figure 5.. BAL Cell Type Differential.

Single cell RNA sequencing data from The Human Lung Cell Atlas was used to calculate cell-type specific marker genes and impute cell-type fractions and enrichment scores (CiberSortX). (A) Highly correlated cell types were identified by calculating Spearman rank-based correlation of cell type fractions followed by hierarchical clustering of correlation distances. (B) Cell type enrichment scores were tested for association with pulmonary restriction (No restriction, n=69; restriction, n=35) using Kruskal-Wallis testing.

Figure 6.. Survival Estimates by Pre-HCT BAL Metatranscriptome and Pre-HCT FVC.

A multivariable Cox regression model for the outcome of post-HCT all-cause mortality was fit for all patients (n=104) using pre-HCT FVC_%pred, pre-HCT BAL RNA mass of Haemophilus (vst transformed), and pre-HCT BAL enrichment score of MSigDB Hallmark Apical Surface Gene Set with covariate adjustment for patient age and sex. Model-based mortality estimates are plotted according to pre-HCT FVC (100% versus 80% versus 60%) and high, medium, and low Haemophilus RNA mass and Apical Surface gene expression. High, medium, and low refer to 75^th, 50^th, and 25^th percentile within the cohort.