A Different Microbiome Gene Repertoire in the Airways of Cystic Fibrosis Patients with Severe Lung Disease

Abstract

In recent years, next-generation sequencing (NGS) was employed to decipher the structure and composition of the microbiota of the airways in cystic fibrosis (CF) patients. However, little is still known about the overall gene functions harbored by the resident microbial populations and which specific genes are associated with various stages of CF lung disease. In the present study, we aimed to identify the microbial gene repertoire of CF microbiota in twelve patients with severe and normal/mild lung disease by performing sputum shotgun metagenome sequencing. The abundance of metabolic pathways encoded by microbes inhabiting CF airways was reconstructed from the metagenome. We identified a set of metabolic pathways differently distributed in patients with different pulmonary function; namely, pathways related to bacterial chemotaxis and flagellar assembly, as well as genes encoding efflux-mediated antibiotic resistance mechanisms and virulence-related genes. The results indicated that the microbiome of CF patients with low pulmonary function is enriched in virulence-related genes and in genes encoding efflux-mediated antibiotic resistance mechanisms. Overall, the microbiome of severely affected adults with CF seems to encode different mechanisms for the facilitation of microbial colonization and persistence in the lung, consistent with the characteristics of multidrug-resistant microbial communities that are commonly observed in patients with severe lung disease.

Keywords: bioinformatics; cystic fibrosis; lung disease; lung microbiome; metabolic pathways; shotgun metagenomics; virulence genes.

Figure 1

Principal component analysis (PCA) based on the top twenty metabolic patterns. Each number corresponds to a pathway reported in the figure legend, whereas each point corresponds to a different patient. Point shape reflects patient groups.

Figure 2

Differences in metabolic and regulatory pathways detected with HMP Unified Metabolic Analysis Network HUMAnN. Colors indicate statistically significant differences (p-value < 0.05) after Kruskal–Wallis one-way analysis of variance. Values indicate mean and one standard deviation (bars). Standardized abundances (x axis) were calculated as: [x − m(x)]/sd(x), where “SD” is the standard deviation and “m” is the mean value.

Figure 3

Distribution of functional categories obtained with an evolutionary genealogy of genes: non-supervised orthologous groups (eggnog) analysis. Box plots for each eggNOG category found for both groups of patients are reported in the left panel, whereas relative abundances for each group are reported in the right panel (bars represent the average value and error bars indicate the standard error for each measure). Significant differences between normal/mild and severe groups are flagged with an asterisk. Boxes denote the interquartile range (IQR) between the 25th and the 75th percentile (first and third quartiles), whereas the inner line represents the median. Whiskers represent the lowest and highest values within 1.5 times IQR from the first and third quartiles. Outliers are reported using white circles.

Figure 4

Virulence factor distribution across patient groups. Only factors displaying a significant (Student’s t-test p-value ≤ 0.05) diverging distribution are reported. The bars represent the average value for each category whereas the error bars indicate the standard error.

Figure 5

Antibiotic resistance gene differentiate patient groups. Bar charts report the percentage values of antibiotic resistance gene detected. Significant differences are reported with one asterisk (p-value < 0.05). P values were obtained through a Student’s t-test on the number of genes detected. Bars were drawn by computing the average percentage value of each resistance category whereas error bars are reported using standard errors.

Figure 6

Relative influence of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (K1 to K5), gene frequencies (evolutionary genealogy of genes: non-supervised orthologous groups (eggnog) database, E1 to E6), antibiotic resistance genes (Resfam database A1 and A2), and virulence factors (microbial database of protein toxins, virulence factors and antibiotic resistance genes for bio-defense applications (MvirDB), V1) for taxonomic structure of microbial lung community evaluated through “booster regression trees” models. The relative influence of different principal component to each bacterial species detected was clustered with the unweighted pair-group method with arithmetic mean (UPGMA) method based on Pearson’s correlation and a tree was reported on the right side of the plot. Boxes report the main group of taxa detected through a cluster analysis.