Functional lower airways genomic profiling of the microbiome to capture active microbial metabolism

Abstract

Background: Microbiome studies of the lower airways based on bacterial 16S rRNA gene sequencing assess microbial community structure but can only infer functional characteristics. Microbial products, such as short-chain fatty acids (SCFAs), in the lower airways have significant impact on the host's immune tone. Thus, functional approaches to the analyses of the microbiome are necessary.

Methods: Here we used upper and lower airway samples from a research bronchoscopy smoker cohort. In addition, we validated our results in an experimental mouse model. We extended our microbiota characterisation beyond 16S rRNA gene sequencing with the use of whole-genome shotgun (WGS) and RNA metatranscriptome sequencing. SCFAs were also measured in lower airway samples and correlated with each of the sequencing datasets. In the mouse model, 16S rRNA gene and RNA metatranscriptome sequencing were performed.

Results: Functional evaluations of the lower airway microbiota using inferred metagenome, WGS and metatranscriptome data were dissimilar. Comparison with measured levels of SCFAs shows that the inferred metagenome from the 16S rRNA gene sequencing data was poorly correlated, while better correlations were noted when SCFA levels were compared with WGS and metatranscriptome data. Modelling lower airway aspiration with oral commensals in a mouse model showed that the metatranscriptome most efficiently captures transient active microbial metabolism, which was overestimated by 16S rRNA gene sequencing.

Conclusions: Functional characterisation of the lower airway microbiota through metatranscriptome data identifies metabolically active organisms capable of producing metabolites with immunomodulatory capacity, such as SCFAs.

Figure 1:. 16S rRNA gene, Whole genome (WGS) and RNA sequencing:

Background (BKG), Upper Airway (UA) and Bronchoalveolar (BAL) samples were collected via bronchoscopy; 16S rRNA gene, Whole genome and RNA sequencing was performed. (A) A heatmap based on Bray-Curtis distance for the 16S rRNA gene sequencing, illustrates the top taxa for all samples. Hierarchical clustering showed two clear clusters, one with BKG samples and BAL samples similar to BKG (Background Predominant Taxa) and another with UA samples and BAL samples similar to UA (Supraglottic Predominant Taxa). (B) Dirichlet Multinomial Modelling (DMM) showed 2 clusters had the best model fit for the 16S rRNA gene sequencing. (C) α Diversity, measured by Shannon Index, showed significant (Wilcoxon) difference between all samples and lowest diversity in UA and among BAL samples that clustered to BAL.16S.SPT by DMM. (D) Beta Diversity, measured by Bray-Curtis, also indicates a significant (PERMANOVA) difference between all samples for 16S rRNA gene sequencing. (E) Bacterial load, measured by ddPCR showed highest levels in UA samples (Kruskal-Wallis). BAL Samples also had higher levels when compared to BKG samples. (F) The inferred metagenome was assessed using PICRUST highlighting several significantly enriched pathways (colored in red). (G) β Diversity for WGS, measured by Bray-Curtis, showed a significant (PERMANOVA) difference between all samples, with UA samples separate from BKG and BAL.BPT samples. Three BAL.SPT samples clustered with UA Samples. (H) β Diversity for RNA, measured by Bray-Curtis, showed a significant (PERMANOVA) difference between all sample types. Two BAL.SPT samples clustered with UA samples. (I) Z Transformed Bray-Curtis Distance between BAL samples and paired UA samples showed clear separation of BAL.16S.BPT and BAL.16S.SPT samples in 16S rRNA gene sequencing. This separation was not as clear in WGS and RNA.

Figure 2:. Taxonomic annotation of all three sequencing data types.

DESEQ2 analysis of taxonomic annotation (at the genus level) between BAL.16S.SPT versus BAL. 16S.BPT samples (FDR <0.05) was performed on 16S rRNA gene sequencing data (A), WGS data (C) and RNA metatranscriptome data (E). Circle size is representative of relative abundance. Gene Set Enrichment Analysis (GSEA) was used to compare the taxonomic signatures identified as distinctly enriched in BAL.16S.SPT vs. BAL. 16S.BPT samples across the different sequencing platforms (B, D, F).

Figure 3:. Functional annotation of all 3 sequencing data types.

(A) Gene Set Enrichment Analysis (GSEA) comparing functional signatures identified across the different sequence data types as distinctly enriched in BAL.16S.SPT vs. BAL. 16S.BPT samples based on KEGG Orthology (KO) annotation (differential enrichment performed based on DESEQ2 analysis). (B) KOs were summarized to associated pathways and differential expression between BAL.16S.SPT and BAL.16S.BPT are displayed as circles for 16S rRNA gene sequencing, diamonds for WGS and squares for RNA. Coloring indicates statistical significance (DESeq2 p<0.05) for each sequence data type and size is relative to the amount of KOs contributing to that pathway. Two pathways highlighted in red include Fatty Acid Biosynthesis, which shows concordance of directionality between the three sequence data types and Fatty Acid Metabolism, which shows discordance.

Figure 4:. Concentrations of Short Chain Fatty Acid (SCFA) in bronchoscopy samples:

A panel of SCFAs were measured and compared (Kruskal-Wallis) in Background (BKG), Upper Airway (UA) and Bronchoalveolar (BAL) samples by GC-MS. SCFA were derived from the linear phase of the standard curve leading to the following cutoffs values (dotted line): (A) 1μM for Acetate, (B) 0.6 μM for Propionate (C) 0.01 μM for Isovalerate and (D) 1μM for Butyrate.

Figure 5:. Diversity correlations with SCFA measurements:

(A) Levels of SCFAs with Acetate, Propionate, Isovalerate and Butyrate were tested (PERMANOVA) against Beta Diversity distribution of data from all three sequencing techniques in BAL samples. Relative abundance of three KOs with direct annotation to measured SCFAs were compared across sample types: K01738 (Acetate), K00925 (Propionate) and K01034 (Butyrate) with (B) 16S rRNA gene sequencing, (C) Whole Genome Sequencing and (D) RNA metatranscriptome sequencing. (E) RNA metatranscriptome taxonomic annotation for these three SCFAs-associated KOs in UA, BAL.RNA.SPT, BAL.RNA.BPT and BKG samples are represented here. Each circle represents a different sample type and colors indicate a different taxonomic annotation.

Figure 6:. Mouse experiment with 16S rRNA gene and RNA metatranscriptome sequencing:

(A) Visual schematic of the experiment, mice (n=17) were inoculated with a mixture of Prevotella, Streptococcus and Veillonella (MOC) and sacrificed at specific time intervals: 1 Hour, 4 Hours, 1 Day, 3 Days, 7 Days. BAL samples were analyzed by 16S rRNA gene sequencing (B-D) and RNA metatranscriptome sequencing (E-G): Principle coordinate analysis was performed with Bray-Curtis Distances by time point (B, E). Mean inter-group distance between sample time point and PBS was calculated (C, F). Relative abundance for taxa annotated to Prevotella, Streptococcus and Veillonella were calculated for each time point (D, G)