ABSTRACT
Rationale Impaired baseline lung function is associated with mortality after allogeneic hematopoietic cell transplantation (HCT). Limited knowledge of the molecular pathways that characterize pre-transplant lung function has hindered the development of lung-targeted interventions.
Objectives To elucidate the biologic and microbiologic correlates of impaired lung function in pediatric allogeneic HCT candidates.
Methods Between 2005-2016, 104 patients with malignant and non-malignant disorders ages 4-19 years underwent paired pulmonary function testing (PFT) and bronchoalveolar lavage (BAL) a median of 1-2 weeks prior to allogeneic transplant in Utrecht, the Netherlands. Cryopreserved BAL underwent RNA sequencing followed by alignment to microbial and human reference genomes for microbiome and gene expression analyses.
Measurements and Main Results 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. BAL microbiome depletion of commensal supraglottic taxa such as Haemophilus and enrichment of nasal and skin taxa such as Staphylococcus were associated with worse measures of lung capacity and gas diffusion. In addition, impaired lung capacity and diffusion were also associated with gene expression signatures of alveolar epithelial proliferation, epithelial-mesenchymal transition, and downregulated immunity, suggesting a post-injury pro-fibrotic response. Detection of microbial depletion and abnormal epithelial gene expression in BAL enhanced the prognostic utility of pre-HCT PFTs for the outcome of post-HCT mortality.
Conclusions These findings suggest a novel and potentially actionable connection between microbiome depletion, alveolar injury, and pulmonary fibrosis in the pathogenesis of pre-HCT lung dysfunction.
INTRODUCTION
Allogeneic hematopoietic cell transplantation (HCT) is a life-saving cellular therapy used to cure underlying hematopoietic defects such as hematologic malignancies, primary immunodeficiencies, bone marrow failure syndromes, and hemoglobinopathies. While more than 7,500 children receive this treatment annually in the United States and Europe, the success of this therapy is limited by a broad range of post-HCT pulmonary injuries that develop in 10-40% of pediatric recipients in response to chemotherapeutic toxicity, infection, and impaired or dysregulated immunity (1–3). Comprehensive assessment of pulmonary health prior to HCT conditioning is critical in order to identify high-risk patients who may benefit from modified conditioning chemotherapy regimens, closer surveillance, and risk-stratified treatments aimed at mitigating the development of irreversible pulmonary toxicity (4–6).
In children, evaluation of pulmonary health prior to HCT is not standardized but may include symptom screening and physical exam, nasopharyngeal swab for common respiratory viruses, pulmonary function testing (PFT) if developmentally able (typically age ≥6 years), and targeted high-resolution chest tomography (HR-CT) and/or bronchoalveolar lavage (BAL) in patients with identified or suspected abnormalities (7, 8). As a sensitive metric of organ function, abnormal PFTs are detected in up to 60% of pediatric HCT candidates predominantly due to restriction and impaired diffusion (9–11) and are associated with the development of post-HCT complications including interstitial pneumonitis, pulmonary graft-versus-host disease (GVHD), and transplant-related mortality (TRM) (2, 3, 12–16).
To improve outcomes, it remains crucial to identify risk-factors for pre-HCT pulmonary dysfunction and elucidate pathobiologic correlates. To date, known risk-factors include exposure to alkylating and other pulmotoxic chemotherapy (e.g. busulfan, bleomycin, cyclophosphamide, carmustine, fludarabine, lomustine), thoracic radiation (particularly in concert with radiomimetics), thoracic surgery, chronic endothelial injury due to hemolysis and/or transfusion dependence, chronic inflammation, and recurrent pulmonary infections (17–19). Most patients with significantly impaired pre-HCT pulmonary function experience a number of these risk-factors, suggesting a multimodal pathogenesis. Despite the need for deeper understanding of the pathogenesis of pre-HCT pulmonary dysfunction, translational studies are limited by scarcity of investigable pulmonary samples from children.
In this study, we aimed to evaluate associations between pre-HCT pulmonary function and measures of pulmonary biology as assessed by metatranscriptomic RNA sequencing of paired BALs. We recently analyzed pre-HCT BAL samples from a cohort of pediatric HCT candidates and found that pre-HCT pulmonary bacterial depletion, viral enrichment, inflammation, and epithelial activation were associated with the development of post-HCT pulmonary toxicity (20). A subset of these patients were developmentally able to perform contemporaneous PFTs, presenting an opportunity to link the domains of organ function, biology, and microbiology within a cohort of high-risk patients. We hypothesized that the combination of pulmonary microbiome data paired with gene expression data could be associated with abnormal PFTs prior to HCT, thus providing further evidence for the pathobiologic role of the pulmonary microenvironment in pediatric HCT candidates.
METHODS
Study Design and Patients
As previously reported, pediatric allogeneic HCT candidates ≤19 years of age with malignant and non-malignant diseases underwent protocolized evaluation of pre-HCT pulmonary health between 2005-2016 in Utrecht, the Netherlands (11). This study included the subset who completed both pre-HCT PFTs and pre-HCT BAL prior to the start of conditioning chemotherapy. Informed consent was obtained with Institutional Review Board approval (#05/143 and #11/063). Details of approach to peri-HCT clinical care are described in Supplemental Methods 1.
PFT Measurements
PFTs included spirometry, whole body plethysmography, and carbon monoxide diffusion capacity and were performed in developmentally capable children ages 4 and older according to established criteria (21, 22). Pre-bronchodilator PFT values of adequate performance quality were normalized for age, height, sex, and race and expressed as (1) a percentage deviation from the mean (% predicted) according to healthy Dutch children and (2) a standardized deviation from the mean (z-score) according to the Global Lung Index (GLI) (23, 24). Restriction, obstruction, mixed restriction/obstruction, diffusion impairment, and air trapping were identified first using conventional PFT cutoffs and were compared to GLI-recommended cutoffs wherein PFTs below the 5th percentile (z-score <-1.64) were considered below the lower limit of normal (Supplemental Methods 2) (21, 24).
BAL Measurements
BAL was performed in the absence of lower respiratory tract symptoms and at the time of another clinically-indicated procedure (e.g. bone marrow aspirate, central line placement) (19, 25). Cryopreserved aliquots of BAL underwent metatranscriptomic RNA sequencing as previously described (20). Briefly, samples underwent mechanical homogenization followed by RNA extraction and purification (26, 27). Sample RNA was combined with control spike-in RNA and underwent reverse transcription, library preparation, and 125nt paired-end sequencing on a NovaSeq 6000 instrument to a target depth of 40 million read-pairs per sample (Supplemental Methods 3). Resultant .fastq files underwent quality filtration followed by parallel alignment to both human and microbial reference genomes using the IDseq pipeline (28). Microbe count data were adjusted for laboratory contamination and estimated microbe mass was calculated using the reference spike-in (Supplemental Methods 4) (29).
Data Analysis
Associations between PFTs and patient characteristics were tested using Kruskal-Wallis tests and Spearman correlations for categorical and continuous traits, respectively. Associations between pre-HCT PFTs and post-HCT outcomes were tested with Kaplan-Meier estimates and Cox proportional hazards models for all-cause mortality and cumulative incidence function with competing risks regression for fatal lung injury. The pulmonary microbiome was described using taxa correlation matrices and principal component analysis (PCA). Associations between PFTs and microbial masses were tested using negative binomial generalized linear models (NB-GLM) in the R statistical environment (edgeR) (30). Gene expression was characterized by calculating enrichment scores for the Hallmark Gene Sets in the Molecular Signatures Database (MSigDB) (31). Associations between PFTs and differentially expressed genes were tested using NB-GLM and characterized using pathway analysis (enrichR) (32). BAL cell type proportions and cell type-specific gene expression profiles were imputed using CiberSortX and the Human Lung Cell Atlas (Supplemental Methods 5) (33, 34).
Data Sharing Statement
Sequencing files are available through dbGaP study accession phs002208.v1.p1.
RESULTS
Patient characteristics
Contemporaneous pre-HCT PFTs and BAL fluid were available for 104 patients. PFTs and BAL were typically performed 1-2 days apart and 1-2 weeks prior to HCT conditioning. Most patients were between 8-19 years of age, of self-reported Northern European ancestry, and preparing to undergo 1st allogeneic HCT for underlying hematologic malignancy (Table 1). Most patients received myeloablative conditioning chemotherapy followed by an umbilical cord blood allograft. At 2-year post-HCT follow-up, 33 of 104 patients (31.7%) had died from relapsed disease or transplant toxicity.
Impaired Pulmonary Function Prior to HCT Is Common
All patients performed spirometry and a subset of typically older children performed plethysmography (n=78) and diffusion testing (n=63, Table 2). To determine the inter-relatedness of the different PFT measurements, we constructed a correlation matrix which showed broad grouping of PFTs into 4 clades pertaining to lung capacity, obstruction, diffusion, and air trapping (Figure E1). Using conventional cutoffs of absolute ratios and percentage deviations from the population mean (% predicted), PFT abnormalities were detected in 80 of 104 patients due to restriction (n=35/104), obstruction (n=12/104), mixed restriction/obstruction (n=1/104), impaired diffusion (n=21/63), and air trapping (n=44/78). Using GLI cutoffs based on standard deviations from the population mean (z-scores), PFT abnormalities were detected in 54 of 104 patients due to restriction (n=29/104), obstruction (n=7/104), impaired diffusion (n=14/63), and air trapping (n=16/78, Table E1). As GLI z-scores are not calculable for all measured PFTs, we used PFTs expressed as percentage of predicted for subsequent analyses.
Abnormal Pulmonary Function Varies According to Age and Sex
We tested for clinical characteristics associated with PFT abnormalities and identified that older age was associated with restrictive lung disease as well as worse measures of obstruction (e.g. MEF25, MEF50, MEF75, and PEF, Figure 1a). Female sex was associated with worse measures of both lung capacity (e.g. FVC, Va) and diffusion (e.g. DLCO/Va) but better measures of obstruction (FEV1/FVC, Figure 1b). We tested for but did not identify an association between HCT indication and pre-HCT PFTs (Supplemental Data).
Abnormal Pulmonary Function Is Associated with HCT Outcomes
Previous studies have shown impaired PFTs to be associated with post-HCT mortality (2, 3, 12–16). We confirmed these findings by demonstrating an association between reduced pre-HCT lung capacity and the development of all-cause mortality after HCT (HR 1.21 per 10% decline in FVC%pred, 95% CI 1.03-1.43, p=0.024, Figure 1c). To test whether these poor outcomes were due to acute lung injury specifically, we calculated cumulative incidence with competing risks regression and found that FVC%pred was associated with the development of fatal lung injury in the first 180 days post-HCT (p<0.001, Figure E2).
BAL Microbiome Composition Varies With Pulmonary Function
Our prior studies have associated alterations in the pre-HCT pulmonary microbiome composition with the later development of post-HCT outcomes (20). To determine whether the composition of the pre-HCT pulmonary microbiome is associated with lung function at the time of BAL measurement, we compared microbiome findings to contemporaneously performed PFTs. We first surveyed the composition of BAL microbiomes and found two oppositely correlated clusters of microbes: one containing common supraglottic taxa (e.g. Haemophilus, Streptococcus, Neisseria) and a second inversely correlated cluster containing common nasal, skin and environmental taxa (e.g. Staphylococcus, Corynebacterium, Figure E3). We represented these clusters mathematically using PCA and found that lower coordinates along principal component 1 (with top contributions from oropharyngeal taxa Prevotella, Veillonella, Rothia, and Streptococcus) were associated with worse measures of lung capacity (e.g. FVC, Figure E4, Supplemental Data).
Next, we used NB-GLM to identify specific taxa associated with each PFT. In models adjusted for patient age and sex, we found that lower BAL quantities of supraglottic taxa (e.g. Haemophilus, Neisseria) and greater BAL quantities of skin/environmental taxa (e.g. Staphylococcus, Saccharomyces, Corynebacterium) were associated with worse measures of lung capacity (FVC), obstruction (MEF25), diffusion impairment (DLCO), and air trapping (RV/TLC) (Figure 2, Figure E5, Supplemental Data). In addition, greater pulmonary fungal mass was associated with worse measures of lung capacity (e.g. FVC) and diffusion (e.g. DLCO), and greater mass of the Actinobacteria phylum, which includes the skin-associated genera Corynebacterium, Micrococcus, and Cutibacterium, was associated with worse measures of lung capacity (e.g. FVC) and obstruction (e.g. MEF25, MEF75). In contrast, lower mass of the Fusobacterium phylum was associated with worse measures of lung capacity (Va) and diffusion (DLCO). We did not find an association between PFTs and total viral RNA mass (Supplemental Data).
Some microbiome studies have suggested that aggregate microbiome traits such as alpha diversity, dominance, and richness, are associated with severity of diseases such as COPD and IPF (35, 36). In this cohort, neither alpha diversity nor richness were associated with PFTs. There was a trend towards an association between greater percentage of the microbiome occupied by the most abundance taxa (i.e. dominance) and lower measures of lung capacity (e.g. FVC, p=0.071).
BAL Gene Expression Is Correlated with Pulmonary Function
To determine whether impaired pulmonary function was associated with unique gene expression signatures, we first utilized an overview approach and calculated BAL expression enrichment scores for each of the n=50 Hallmark Gene Sets in the Molecular Signatures Database. We found that measures of lung capacity and diffusion were positively associated with enrichment of pathways relating to immune activation, metabolism, and cell division and were negatively associated with pathways relating to epithelial mesenchymal transition and cell fate/differentiation (Figure 3, Supplemental Data).
We next tested for specific genes whose expression levels were associated with each PFT. Pathway analysis showed that lower levels of both FVC and DLCO/Va were associated with greater expression of pulmonary epithelial genes, including surfactant genes (e.g. SFTPA1/2, SFTPC, SFTPD), keratinization genes (e.g. KRTAP3.2, KRTAP5.6), and epithelial sensory/stimulus response genes (e.g. OR2A1, OR4E2), as well as lower expression of innate immunity genes (e.g. CXCL8, HLA-DRA, IL1B, LYZ; Figure 4a, 4b). In contrast, lower FEV1/FVC was associated with greater expression of innate immunity genes and lower expression of alveolar epithelial genes. These findings persisted after adjusting models for age and sex. Few genes (n=9) were associated with a metric of air trapping (RV/TLC, Supplemental Data).
BAL Cell Fraction is Associated with Pulmonary Function
Since BAL samples are heterogeneous cell mixtures, we posited that these results could be due to either different proportions of cell types within samples, different expression patterns of each cell type within samples, or both. Therefore, we used in silico cell deconvolution techniques based on the Human Lung Cell Atlas to impute lung cell proportions in each sample. BAL cell types clustered into 3 groups consisting of innate immune cells, upper airway cells, and lower airway and adaptive immune cells (Figure 5a). Patients with restriction had greater BAL enrichment scores for type 1 and 2 alveolar epithelial cells (AECs) (Figure 5b). In addition, patients with impaired diffusion had greater BAL enrichment scores for type 2 AECs (Supplemental Data). Next, using the CiberSortX algorithm (33), we estimated the average gene expression profile for each BAL cell type. In patients with restriction, type 1 AECs showed increased transcription of genes related to epithelial-mesenchymal transition, hedgehog signaling, and KRAS signaling, and T- and B-lymphocytes showed an overall decrease in transcription of genes relating to growth and inflammation (Supplemental Data). We did not detect differences in innate immune cell or upper airway epithelial cell-specific gene expression among patients with vs. without restriction.
Interactions Between PFTs, BAL Metatranscriptome, and Post-HCT Mortality
Given prior findings that pre-HCT BAL metatranscriptomes are associated with post-HCT mortality (20), we next sought to determine whether these associations are independent of measured pulmonary function. Using multivariable modeling, we found that lower BAL mass of oropharyngeal taxa (e.g. Haemophilus, Neisseria) was associated with post-HCT all-cause mortality independent of FVC, age, and sex ((Supplemental Data). We identified a statistical interaction whereby progressively lower BAL mass of oropharyngeal taxa strengthened the association between poor pre-HCT FVC and post-HCT mortality. In addition, lower BAL expression of genes related to airway epithelial cell integrity, polarity, and junction formation (e.g. Hallmark Apical Surface) were also associated with post-HCT mortality independent of FVC, age, and sex. The addition of BAL Haemophilus mass and BAL Hallmark Apical Surface gene set enrichment to a Cox survival model of FVC, age, and sex significantly improved model performance across strata of lung function (likelihood ratio test p=0.002, Figure 6).
DISCUSSION
The primary finding of this study is that in pediatric HCT candidates, pulmonary microbiome and gene expression signatures are associated with pulmonary function and may augment the utility of PFTs in the stratification of post-HCT mortality risk. Specifically, we found that the pulmonary microbiomes of patients with restriction and diffusion impairment were depleted of common supraglottic taxa and instead enriched in skin, nasal, and environmental taxa. In addition, the pulmonary transcriptomes of patients with restriction and diffusion impairment showed increased activation of alveolar epithelial cells and a deficit in immune cell signaling consistent with post-injury fibrotic transformation. Further, BAL signatures of commensal microbial depletion and diminished epithelial homeostasis were associated with post-HCT mortality independent of pre-HCT pulmonary function. Together these findings suggest a novel connection between lung function, biology, and microbiology in pediatric HCT candidates and suggest possible biologic targets to optimize lung function in children preparing to undergo HCT.
First, our finding that pulmonary microbiome depletion is associated with impaired lung function in pediatric HCT candidates is novel and supports our previous finding that pre-HCT pulmonary microbiome depletion is associated with post-HCT lung injury (20). Although studies of the pulmonary microbiome in children are rare, investigations of the more readily-studied intestinal microbiome have similarly shown that loss of intestinal biodiversity and biomass are associated with organ dysfunction as well as poor post-HCT outcomes (37). Whether pulmonary microbial dysbiosis can directly contribute to lung disease or is simply a biomarker of other exposures remains unclear. On the one hand, microbial dysbiosis is linked to antimicrobial exposure, which is more common in patients with underlying diseases such as malignancy; such patients have confounding reasons for poor lung function such as greater exposure to pulmotoxic chemotherapy. However, emerging studies suggest that microbial dysbiosis disrupts the control of pathogen overgrowth and contributes to a loss of immune tolerance at the epithelial-parenchymal boundary, and therefore reciprocal causality between lung dysfunction and microbiome disturbances seems likely (35, 38). To that end, we found that pulmonary biomass of supraglottic commensals was inversely related to biomass of skin, nasal, and environmental taxa, supporting the theory that maintenance of the respiratory microbiomes may exclude invasion by outside organisms. Ultimately, strategies to modulate the lung microbiome in order to improve pulmonary health remain in development and include investigation of both ingested and inhaled probiotics and other microbial metabolites (39, 40).
A second finding of this work is the association between altered alveolar epithelial gene expression, restriction, and impaired diffusion. By combining differential gene expression with cell type imputation, we identified a greater proportion of type 1 and 2 AECs in patients with restriction and diffusion and cell-type specific gene expression suggested an increased transcriptional state of these cells, associated with growth factors and increased kinase expression. The increased expression of surfactant genes and overall AEC activity is consistent with response to alveolar injury and has been demonstrated in numerous cohorts of lung-injured patients (41, 42). By combining hallmark gene set enrichment with these findings, we found that pathways related to the lung stroma were significantly upregulated, including angiogenesis, myogenesis, epithelial mesenchymal transition, and cell fate (Wnt/β-catenin signaling and hedgehog signaling). Activation of stromal pathways in the lung has been associated with a pro-fibrotic signaling state and has been identified in cohorts of patients with IPF, thus providing a plausible mechanistic link between these gene expression findings and the observed impairment in PFTs (43, 44). As lung-resident mesenchymal stem cells undergoing myofibroblastic differentiation strongly express the Wnt/β-catenin and hedgehog signaling pathways, anti-fibrotic agents targeting this pathway such as pirfenidone remain of utmost interest (45, 46).
A third finding of this work is that patients with pulmonary restriction and impaired diffusion displayed an overall deficit in lower airway immune signaling. In addition to lower overall expression of immune genes among patents with restriction, cell type deconvolution inferred a less active transcriptional state for lymphocytes specifically. It seems counterintuitive that patients with poor PFTs would have a relative deficit in pulmonary immune activity, especially given our prior findings that pre-HCT pulmonary inflammation was associated with post-HCT lung injury. There are three potential explanations. First, the latter finding of pulmonary inflammation was tied to pre-HCT viral infection, which was more prevalent in children <4 years old. As younger children are not typically able to perform PFTs, they were excluded from this present study. Second, while pulmonary inflammation is a hallmark of acute lung injury including chemoradiation-induced lung injury, it typically subsides after the initial insult and is followed by a profound anti-inflammatory/pro-fibrotic phase (47), which more likely represents the phase of injury for the patients in this study. Third, patients with impaired lower airway immunity may be at greater risk for infections and associated sequelae such as post-infectious fibrosis. Ultimately, the significance of impaired pulmonary immune environment in pediatric HCT candidates with poor lung function needs to be explored with novel approaches in the future.
A fourth and unanticipated finding of this work was the identification of significant sex differences in pulmonary function. Although both male and female study participants had worse lung function than age- and sex-matched controls, females had profoundly worse measures of both restriction and impaired diffusion and males had slightly worse measures of obstruction. Long-term follow-up studies of pediatric oncology and HCT patients have also demonstrated associations between female sex, restriction, and impaired diffusion as well as between male sex and obstruction; this study adds to the field by identifying these disparities prior to HCT. A likely culprit for these toxicities are underlying differences in chemotherapy pharmacokinetics and pharmacodynamics, which merit exploration in future studies (48–50).
Our study has several strengths. First, the analysis of contemporaneously measured PFTs and BAL is unique among pediatric cohorts and allows novel insights into lung function, biology, and microbiology. Second, PFT measurements included GLI Z-scores according to updated ATS/ERS recommendations. Third, we analyzed BAL using a paired assessment of microbiome and gene expression. Several limitation merit comment. First, as with all observational human studies, we are unable to prove mechanism or causation among the factors measured and described herein. Second, we were unable to ascertain all possible risk-factors that might have influenced the pulmonary microenvironment, including history of pulmonary infections, prior antimicrobial exposure, and detailed chemotherapeutic regimens. Finally, younger children were necessarily excluded due to inability to perform PFTs. However, given our findings associating PFTs with lower respiratory biology, this study raises the possibility that BAL RNA sequencing could be a useful surrogate to detect pro-fibrotic states in those unable to perform PFTs.
In conclusion, in this study of pediatric HCT candidates, predominant patterns of pulmonary restriction and impaired gas diffusion were associated with microbiome depletion, proliferation of the alveolar epithelium, and a deficit in immune signaling. Both pre-HCT PFTs and BAL metatranscriptomes were independently associated with post-HCT mortality. These findings significantly increase our knowledge of the pathobiology of lung dysfunction in patients preparing to undergo HCT and offer numerous potential biologic targets that might be manipulated to optimize lung health in pediatric HCT patients.
Data Availability
Data Sharing Statement: Sequencing files are available through dbGaP study accession phs002208.v1.p1.
ACKNOWLEDGEMENTS
We thank Dr. Dean Sheppard, MD for his critical review of this manuscript.
Footnotes
↵* Zinter MS and Versluys BA share first authorship.
↵# DeRisi JL and Boelens JJ share senior authorship.
Funding: National Marrow Donor Program Amy Strelzer Manasevit Grant (MSZ), NICHD K12HD000850 (MSZ), NHLBI K23HL146936 (MSZ), American Thoracic Society Foundation Grant (MSZ), NHLBI R01HL134828 (MAM), NHLBI R35HL140026 (MAM), Chan Zuckerberg Biohub (JLD).