Transitions in lung microbiota landscape associate with distinct patterns of pneumonia progression

Abstract

The precise microbial determinants driving clinical outcomes in severe pneumonia are unknown. Competing ecological forces produce dynamic microbiota states in health and disease, and a more thorough understanding of these states has the potential to improve pneumonia therapy. Here, we leverage a large collection of bronchoscopic samples from patients with suspected pneumonia to determine lung microbial ecosystem dynamics throughout the course of pneumonia. We combine 16S rRNA gene, metagenomic, and metatranscriptomic sequencing with bacterial-load quantification to reveal clinically relevant drivers of pneumonia progression. Microbiota states are predictive of pneumonia subtypes and exhibit differential stability and pneumonia therapy response. Disruptive forces, such as aspiration, are associated with cohesive changes in gene expression and microbial community structure. In summary, we show that host and microbiota landscapes change in unison with clinical phenotypes and that microbiota state dynamics reflect pneumonia progression. We suggest that distinct pathways of lung microbial community succession mediate pneumonia progression.

Keywords: infection; metagenomics; metatranscriptomics; multiomics; pneumonia; respiratory tract microbiome.

Figure 1.. Multiomics of the lung microbial ecosystem during pneumonia reveals diverse associations with clinical features.

(A) Clinical characteristics of the SCRIPT cohort based upon clinical outcome. (“−” indicates patients with negative binned clinical outcomes [e.g., patient expires], “+” indicates patients with positive binned clinical outcomes [e.g., patient discharged and sent home]) (B) UpSet plot of multiomics sampling at the same time-point. Colors distinguish sample as either a baseline or follow-up BAL. (AMP = 16S rRNA gene amplicon, MGX = metagenomic, MTX = metatranscriptomic [including host-transcriptomics]) (C) Principle coordinate analysis of the weighted UniFrac distance metric derived from amplicon profiles (genus-level). Colors are indicative of pneumonia category. Gray dots in the background are the shadows of points from all four plots combined. Percentages on axes are the variance explained by the given principle coordinate. See Figure S2A for remaining profiles. (HAP = hospital acquired pneumonia, VAP = ventilator-associated pneumonia, CAP = community acquired pneumonia, NP = critically-ill non-pneumonia control.) (D) Permutational multivariate analysis of variance analysis (PERMANOVA) quantifies the amount of variance in distance space explained by a given metadata features (e.g., pneumonia diagnosis) and tests for significance association. Percentages and color represent variance explained (R²). Columns are the different multiomic profiles. Bracketed numbers on right of y-axis metadata labels represent degrees of freedom; numerical variables have a degree of freedom of one. Significant association with high variance explained indicates metadata features as drivers of variation in the gene expression or microbiota landscape. 16S rRNA gene amplicon sequencing profiles include: genus-level abundances (AMP [Genus]) and ASV-level abundances (AMP [ASV]); shotgun metagenomic profiles include: microbial gene family (DNA [KO]), species-level (DNA [Taxonomy]), and viral abudances (DNA [vOTU]); metatranscriptomic profiles include: alveolar macrophage gene expression (RNA [Host Transcriptomics]), genus-level microbial expression (RNA [Taxonomy]), and microbial gene-family expression (RNA [KO]); and respiratory culture profiles include: genus-level cultivation (CFU [Culture], presence/absence). (* FDR P < 0.05, ** FDR P < 0.01, *** FDR P < 0.001; PEEP = positive-end expiratory pressure, FiO2 = fraction of inspired oxygen, Binned Clinical Outcome = positive or negative discharge status as in [A]). (E) Shannon diversity of different multiomics profiles and genus-level culture data. (Boxplot configuration: Center line = median, box limits = upper and lower quartiles, whiskers = 1.5x interquartile range, points = outliers.)

Figure 2.. Pneumonia infection associates with an altered microbiota landscape indicative of aspiration-mediated disruption.

(A) Abundance of the top differentially abundant (FDR P < 0.05) genes, species, and bacterial genera in each pneumonia category (i.e., HAP, VAP, CAP) relative to NP. Bar plots are the proportion of samples with zero-count therefore showcasing feature prevalence; bars are scaled such that touching the correspondingly colored line above indicates the feature was undetected in all samples for that group. Kernel distributions were calculated based on the subset of samples with detectable abundance after centering by the median and log₂ transformation; heights are scaled by the proportion of detectable samples. Genes include their KEGG orthology term. (* FDR P < 0.05, ** FDR P < 0.01, *** FDR P < 0.001). (B) Distribution of the mean dissimilarity to non-pneumonia (MDNP) score quantifying microbiome disruption relative to non-pneumonia control group. Score is calculated using the weighted UniFrac distance from amplicon profiles. The 90^th percentile of MDNP score within NP was used as a threshold to determine microbiota disruption in patients with pneumonia. (C) Abundance of the top differentially abundant (FDR P < 0.05) genes, species, and bacterial genera in disturbed microbial communities (>90^th) relative to communities with structure typical of NP (<90^th). Microbiome samples were binned into typical and disturbed subsets based on the 90^th percentile of MDNP score. Above this threshold, there is a 10% chance of a patient without pneumonia to have that particular arrangement of microbiota. (D) Relationship between bacterial load, amylase activity, and MDNP score. Shaded region represents 95% confidence interval. Statistics show Spearman rank correlation test. (E) Top differentially abundant genera (amplicon); samples ordered by increasing levels of amylase activity. (F) Distribution of bacterial load, amylase activity, and MDNP score binned by culture results and antibiotic usage at time of BAL. Stars represent statistical significance as determined by Wilcoxon test. (* FDR P < 0.05, ** FDR P < 0.01, *** FDR P < 0.001; Boxplot configuration: Center line = median, box limits = upper and lower quartiles, whiskers = 1.5x interquartile range, points = outliers.). Acronyms: HAP = hospital acquired pneumonia, VAP = ventilator-associated pneumonia, CAP = community acquired pneumonia, NP = critically-ill non-pneumonia control.

Figure 3.. An inflammatory response gene expression pattern in alveolar macrophages is elevated with increasing microbial community disruption.

(A) Volcano plot showing results of host differential expression analysis using DESeq2 (N total = 106 samples). Red points indicate genes with FDR P < 0.05 and log₂ fold change with an absolute value greater than 0.5 (n=66 differentially expressed host genes). Full table of results available in Table S4. (B) Relationship between NINJ, PFKFB3, IL1B, and NFKB2 host gene expression with MDNP score (genes also in C). Shaded region represents 95% confidence interval. (C) Heatmap showcasing the top 20 most significant differentially expressed host genes. Annotation bar is MDNP score; columns in heatmap distinguish samples that are greater than the 90th percentile of MDNP score. Heatmap values are the log₂ of normalized expression (total sum scale of RPKM counts with pseudocount of 1). (D) Gene ontology enrichment analysis of host biological processes upregulated with increasing MDNP score. Top 15 enriched pathways selected by highest fold enrichment (FDR P < 0.05, count of at least 5 differentially expressed host genes in that pathway). Full table of GO term enrichment results available in Table S4.

Figure 4.. A posteriori identification of pneumotypes suggests stabilizing selective forces canalize community structure.

(A) Ordination of weighted UniFrac distance derived from genus-level amplicon profiles. Colors represent the different microbiota states of the distal lung (i.e., pneumotypes) identified using cluster analysis. Percentages represent variance explained by the given principle coordinate axis. (B) Summary heatmap visualizing the mean cluster consensus score. Consensus clustering implementing the partition around medoids cluster algorithm was performed to determine number of groups evident in the weighted UniFrac distance space. (C) Trade-offs in diversity (Wilcoxon test) and (D) core phyla differentiate pneumotypes. (E) Most abundant taxa distinguish microbiota states. Taxa with a mean normalized abundance greater than 0.05 were selected (n=12). Stars represent the adjusted p-value of differential abundance analysis comparing pneumotype_M, pneumotype_SP, and pneumotype_OL to pneumotype_SL. (F) Bacterial biomass, (G) amylase levels, (H) MDNP score, (I) and neutrophil abundance differ significantly between microbiota states. Pairwise comparisons show results of Wilcoxon test with Benjamini-Hochberg correction. (J) Heatmaps visualizing pneumotype frequency across pneumonia category (limited to baseline BAL) and (K) clinical outcome (includes baseline and follow-up BAL). Numbers in heatmaps are the count of BAL in each section; color of tiles is the proportion for that column. Stars represent the adjusted p-value of two-sided pairwise exact binomial tests used to determine deviations from expected distributions (i.e., evenly distributed across the column). Pneumonia therapy outcome is categorized as successful (+), indeterminate (+/−), and unsuccessful (−) treatment. (L) Frequency of transitions between pneumotypes. Nodes (circles) represent the different pneumotypes and the circle size is scaled to the number of samples. Edges (arrows) represent transitions between pneumotypes. Edge labels are the frequency of transitions between pneumotypes accounting for transition to outcome (i.e., final BAL are counted as transitioning to clinical outcome rather than to any pneumotype). (* FDR P < 0.05, ** FDR P < 0.01, *** FDR P < 0.001; Pneumonia diagnosis: HAP = hospital acquired pneumonia, VAP = ventilator-associated pneumonia, CAP = community acquired pneumonia, NP = critically-ill non-pneumonia control; Pneumotypes: SL = skin-like, M = mixed, SP = Staphylococcus predominant, OL = oral-like; Pneumonia outcome: − = unsuccessful treatment, + = successful treatment, +/− = indeterminate treatment; Boxplot configuration: Center line = median, box limits = upper and lower quartiles, whiskers = 1.5x interquartile range, points = outliers.)

Figure 5.. Alveolar macrophage gene expression is partially distinct between skin-like and oral-like pneumotypes.

(A) Variation in the alveolar macrophage transcriptome associates with pneumotype_OL and pneumotype_SL. Permutational multivariate analysis of variance analysis (PERMANOVA) performed by comparing host gene expression profiles of a given pneumotype to the remaining samples with known pneumotypes (N total = 106 samples, N pneumotype_SL = 45 samples, N pneumotype_M = 25 samples, N PneumotypSP = 15 samples, N pneumotype_OL = 19 samples). Percentages and color represent variance explained (R2). (* FDR P < 0.05, ** FDR P < 0.01, *** FDR P < 0.001) (B) Volcano plots showing results of host differential gene expression analysis comparing pneumotype_SL (top) and pneumotype_OL (bottom) to remaining samples with known pneumotypes. Analysis performed using DESeq2 with raw read counts. Red points indicate genes with FDR P < 0.05 and log₂ fold change with an absolute value greater than 0.5 (n=99 host genes overexpressed in pneumotype_OL and n=47 host genes underexpressed in pneumotype_SL). Full table of results available in Table S4. (C) Venn diagram showing overlap fraction of host genes underexpressed in pneumotype_SL and overexpressed in pneumotype_SL from figure 5B. (D) Zero log ridge plot showing the top 10 most differentially expressed host genes by FDR P in pneumotype_SL and pneumotype_OL, respectively. (E) Gene ontology enrichment analysis of host biological processes underexpressed in pneumotype_SL (top) and overexpressed in pneumotype_OL (bottom). Top 15 enriched pathways selected by highest fold enrichment (FDR P < 0.05, count of at least 5 differentially expressed host genes in that pathway). Full table of GO term enrichment results available in Table S4.

Figure 6.. The lung microbial ecosystem is complex and rich with interactions across ecosystem components.

(A) Overview of data types and their appropriate correlation metrics used to generate multiomics integration network with HAllA. Multiomics data integration includes feature profiles from four data types: shotgun metagenomic (taxonomic, viral, microbial gene family abundance profiles), 16S rRNA gene sequencing (ASV-level abundance profiles), and metatranscriptomics (host gene, microbial gene family, and microbial genera expression profiles). Clinical variables, binary culture results, and binary antibiotic administration are also included in the analysis. (B) Number of clusters (left) and number of significant feature association pairs (right) for each omics comparison with HAllA. (C) Network visualization of associated omics features identified from HAllA. Edges are associations colored by Spearman rank correlation or normalized mutual information if connecting to a clinical variable node (parallelogram). Nodes are data features (e.g., genera, genes). Node shape represent different profiles. Dotted circles annotate approximate network neighborhoods (3 largest) identified using GLay community structure analysis; colors selected for aesthetic contrasts and circle sizes are not meaningful. Top significant associations from each dataset comparison are visualized (FDR P < 0.05). Nodes with at least 25 significant associations are larger. Nodes high in one pneumotype are colored accordingly. Features that were high in multiple groups or no groups were kept as gray. Full network association information available in Table S5. Full results of community analysis of largest visualized network component available in Table S5. See Figure S4E for a fully labeled network diagram. (D) Selected relationships highlighted from cyan-circled community, (E) from yellow-circled community, (F) and from magenta-circled community (C). Shaded region represents 95% confidence interval.