A novel approach for combining the metagenome, metaresistome, metareplicome and causal inference to determine the microbes and their antibiotic resistance gene repertoire that contribute to dysbiosis

Abstract

The use of whole metagenomic data to infer the relative abundance of all its microbes is well established. The same data can be used to determine the replication rate of all eubacterial taxa with circular chromosomes. Despite their availability, the replication rate profiles (metareplicome) have not been fully exploited in microbiome analyses. Another relatively new approach is the application of causal inferencing to analyse microbiome data that goes beyond correlational studies. A novel scalable pipeline called MeRRCI (Metagenome, metaResistome, and metaReplicome for Causal Inferencing) was developed. MeRRCI combines efficient computation of the metagenome (bacterial relative abundance), metaresistome (antimicrobial gene abundance) and metareplicome (replication rates), and integrates environmental variables (metadata) for causality analysis using Bayesian networks. MeRRCI was applied to an infant gut microbiome data set to investigate the microbial community's response to antibiotics. Our analysis suggests that the current treatment stratagem contributes to preterm infant gut dysbiosis, allowing a proliferation of pathobionts. The study highlights the specific antibacterial resistance genes that may contribute to exponential cell division in the presence of antibiotics for various pathogens, namely Klebsiella pneumoniae, Citrobacter freundii, Staphylococcus epidermidis, Veilonella parvula and Clostridium perfringens. These organisms often contribute to the harmful long-term sequelae seen in these young infants.

Keywords: antibacterial resistance; causal Bayesian network; multi-omics; origin of replication; peak-to-trough ratio (PTR).

Fig. 1.

MeRRCI – an approach to infer consequences of antibiotic resistance from whole metagenomic studies. (a) A high-performance computing framework for computing microbiome-wide eubacterial metagenomic composition (EMC), eubacterial replication rates as the peak to trough ratio (PTR), and profiles of antibiotic resistance (ABR) genes. The reference collection of genome sequences is partitioned into smaller subcollections. Mapping reads to the database of genome sequences and resistance genes are achieved by using bowtie2, resulting in alignment files in SAM format. The setup is well suited for parallelization and the use of high-performance computing environments. The output of this phase includes the three omics profiles – EMC, PTR and ABR. (b) The tool bnlearn can infer causal relationships between the variables of interest. The inputs to this tool are the omics profiles from the previous phase and a set of constraints with the list of allowable edges. These constraints incorporate our biological understanding of the microbiome and its interactions with the host. The output of this phase is a causal graph where the directed edges are causal relationships that were inferred computationally. The constraints are pictured as a network model containing edges that are allowed to be present. Note that different shapes are used for different types of nodes. The allowed types of edges are shown numbered 1–5. Edges of type 1 connect nodes representing taxon relative abundance; edges of type 2 go from replication rates (PTR) of taxa to their relative abundance; edges of type 3 lead from gene expression to replication rate nodes; type 4 nodes connect gene abundance and taxa; type 5 edges connect miscellaneous variables to all other types of nodes. All other edges are disallowed and included in a blacklist.

Fig. 2.

Microbial metagenomic composition (EMC) and replication rate (PTR) profiles. (a) Average re-normalized relative abundance for the most frequently present species found in control (CC) and antibiotic (AC) cohorts. (b) Average replication rates (PTR) values for bacterial species. Statistical significance was evaluated with a Mann–Whitney U test, where: *P≤0.05 and **P≤0.01. Black lines on each bar indicate standard deviation. (c) Spearman correlation coefficients between PTR and relative abundance values of the same taxon computed separately for the two cohorts, CC and AC. The colour of the filled circles signifies the sign of the correlation (blue for positive and red for negative). In contrast, the size of the filled circles represents the absolute value of the correlation coefficient. Crossed entries represent correlations that are not significant (P>0.05).

Fig. 3.

ABR resistance profile. (a) Venn diagram of the number of genes associated with each of the two cohorts (CC and AC) with P≤0.05. (b) Reads per kilobase million (RPKM) for the most abundant ABR gene families found in the dataset, where *P≤0.05 and **P≤0.01. Black lines on each bar indicate standard deviation. Statistical significance for (a-b) were established with Mann–Whitney U tests.

Fig. 4.

Distribution of PTR values. Violin plots showing the distribution of PTR values for taxa in the samples from the control cohort (CC) infants and infants who recently received antibiotic treatment (AC). Mean PTR in each treatment group was compared to CC samples using the Mann–Whitney U test (P-values are as shown in the figure).

Fig. 5.

PTR and microbial relative abundance. Box-and-whisker plots of values of bacterial taxa that significantly changed their replication rates after exposure to TIM or AMP/MEM antibiotics. Significant pairs identified with the Mann–Whitney U test are marked with asterisks (*P≤0.05, **P≤0.01).

Fig. 6.

ABR genes associated with TIM exposure. The ABR genes present in the samples from infants treated with TIM. Each circle includes the total number of ABR genes identified at each analysis step, while the text above the blue arrows summarizes the analysis performed at those steps: 1) PLS-DA to identify genes that are differentially abundant between AMP/MEM vs. TIM cohorts; 2) Mann-Whitney U tests and Benjamini-Hochberg P-value correction for each gene; 3) Selecting genes positively correlated with TIM; 4) Selecting genes associated with K. pneumoniae.

Fig. 7.

Causal network for the infant gut microbiome. Node sizes are proportional to the average value. Nodes correspond to eubacterial relative abundance (circles), PTR values (rectangles), the abundance of ABR genes (triangles), antibiotic dosage administered to infants (pink pentagons) or mothers (yellow pentagons), and various environmental factors (diamonds).

Fig. 7.

Causal network for the infant gut microbiome. Node sizes are proportional to the average value. Nodes correspond to eubacterial relative abundance (circles), PTR values (rectangles), the abundance of ABR genes (triangles), antibiotic dosage administered to infants (pink pentagons) or mothers (yellow pentagons), and various environmental factors (diamonds).

Fig. 7.

Causal network for the infant gut microbiome. Node sizes are proportional to the average value. Nodes correspond to eubacterial relative abundance (circles), PTR values (rectangles), the abundance of ABR genes (triangles), antibiotic dosage administered to infants (pink pentagons) or mothers (yellow pentagons), and various environmental factors (diamonds).