Partial recovery of microbiomes after antibiotic treatment

ABSTRACT

Antibiotics profoundly affect the gut microbiome and modulate microbial communities. We recently observed that antimicrobial drugs also impact the abundance and distribution of antibiotic resistance genes. In this addendum, we reanalyze our ∼1 trillion nucleotide shotgun metagenomic dataset to quantify comprehensive genomic differences at the sequence level before and after antibiotic treatment. We show that 7 day exposure to cefprozil leads to a statistically significant loss of metagenome sequences. Recovery of gut microbiomes 3 months after antibiotherapy was characterized by the emergence of new genome sequences not observed prior to antibiotic exposure. Participants with low initial gut microbiome diversity had an increased amount of sequences related to antibiotic resistance. Therefore, we suggest that while the taxonomical composition of microbiomes is partially affected by the antibiotic, the genomic content and population structure of bacterial communities is noticeably impacted.

KEYWORDS:

• antibiotic resistance
• bioinformatics
• metagenome comparison
• metagenomics
• Microbiome

Introduction

Metagenomic studies demonstrated an uneven recovery of the human gut microbiome after treatment with antibiotics.^1-3 Researchers have even suggested that early exposure to antibiotics can impact future health.^4-6 In our recent study, we explored how a widely prescribed antibiotic, used to cure mild infections, can alter the human gut microbiome.⁷ Our goal was to evaluate the impact of antibiotic treatment on the gut microbiome. To test this, we gave a 7 day course of cefprozil to 18 healthy volunteers and sequenced the microbial content of stool samples before, at the end of the 7 days course of antibiotic treatment and 90 days later. Six control volunteers, receiving no antibiotic, were also enrolled and analyzed at the same time points (day 0, day 7 and day 90). We then evaluated the impact of cefprozil on the gut microbiome of healthy individuals and derived 3 conclusions. By comparing the taxonomical content of stool samples through time and between individuals, we observed that a limited number of bacteria were consistently affected in all treated participants. Only three genera were significantly increased while a dozen other genera were decreased by antibiotic exposure (Fig. 1). In fact, all the taxa significantly and uniformly affected by the antibiotic had an average abundance below 10%. Therefore, our first conclusion was that the antibiotic had a reproducible effect mostly visible on the lower abundance taxa. Our second conclusion was that the initial composition of the gut microbiome was linked to the effect of the antibiotic on this microbiome. Indeed, exposed participants with lower initial microbiome diversity at day 0 had blooms of Enterobacter cloacae complex bacteria at day 7. Our third conclusion was that the antibiotic treatment modulated the gut microbiota in a fashion that led to the apparition of resistance genes that were not detected prior to treatment.

Figure 1. Modulation of bacterial genera by cefprozil. Average change in abundance ranking between day 7 and day 0 samples for exposed participant. Only the 13 genera for which relative abundance was statistically significant between day 7 and day 0 after antibiotics treatment are displayed (p < 0.05 after Benjamini & Yekutieli false discovery rate correction, see Raymond et al.⁷ for details). Error bars represent confidence interval at 95%.

The initial analysis of our shotgun metagenomics data set (∼1 trillion nucleotides) was performed using taxonomic profiling based on k-mers using the Ray Meta software.⁸ In our original publication, we used taxonomical profiling and gene annotation to provide critical information on the impact of antibiotics on the gut microbiome. The results were interpreted with a microbiology perspective, allowing the generation of new hypotheses for better management of antibiotics. This methodology had 2 main limitations. First, it reduced microbial communities to categories of microbes and did not differentiate bacteria of the same taxa that may harbor wholly different sets of genes.⁹ Second, it is biased by the genome sequences available at the time of analysis and therefore will not consider genomes absent from the reference database.

In this paper, we focused on how microbiomes were affected as a whole and performed quantitative analysis of DNA sequence content before and after antibiotic treatment. Indeed, we generated a taxonomy-independent quantitative measure of the similarity between microbiomes to evaluate the impact of exposure to antibiotics. To do so, we compared the k-mer content of assembled metagenomes using our new software implementation, Ray Surveyor (https://github.com/zorino/ray).

The Ray Surveyor software was built to easily perform large comparative genomics analyses using k-mers as the comparison unit. Briefly, it compares the k-mer content of genome assemblies to measure the similarity between samples at the DNA sequence level. K-mers are substrings of DNA of length k. We use a k-mer length of 31 nucleotides because it is sufficiently specific to allow confidence in genomic comparisons.⁸ Ray Surveyor identifies k-mers shared between 2 or more samples by calculating the presence/absence of k-mers in all the samples. The distinctive function of Ray Surveyor is the capability of comparing genome sequences either based on their complete k-mer content or on a specified subset of k-mers, such as an annotation dataset.

Comparison of microbiomes using k-mers

We therefore compared assembled metagenomes (Accession PRJEB8094 in the European Nucleotide Archive) using Ray Surveyor in order to precisely quantify the short-term and long-term impact of a 7 day course of cefprozil on the gut microbiome of each individual without a priori and biases introduced by taxonomical classification. Overall, the genomic content of antibiotic-exposed microbiome assemblies was more different from their initial state than the content of controls without antibiotics (Fig. 2A). Indeed, we observed an average concordance (Jaccard similarity index) of 57.7% in average between day 7 and day 0 for control patients versus 38.4% similarity for patients that were exposed to the antibiotic cefprozil (p = 3.5e-7, Welch Two Sample t-test). These results indicate that the effect of the antibiotic is greater than the temporal metagenomic variability of controls during the course of a week. The measure of unique k-mers lost or gained after antibiotic treatment confirms a loss of genomic diversity after the use of antibiotics compared to controls (Fig. 3A). We then quantitatively assessed how much microbiomes reverted back to their initial sequence composition after 3 months. The Jaccard similarity index between day 90 and day 0 was 56.0% for controls, which was similar to the difference between day 7 and day 0 (Fig. 2B). For exposed participants, the similarity of 43.8% between day 90 and day 0 represented a significant decrease (p = 0.001, Welch Two Sample t-test) but indicated a partial reversion to the initial metagenome composition. In most samples, there is still a marked difference between the controls and the exposed participants at day 90, suggesting that gut microbiomes are more different from their original sequence composition when they were subjected to cefprozil. However, we observed that this difference in exposed participants is now driven by the emergence of sequences not observed at day 0 rather than by loss of k-mers at day 90 (Fig. 3B). This new analysis sketches a complementary yet congruent portrait of our original analysis, which was based on taxonomical profiling. Indeed, in our initial study, we observed that most of the participants returned to their initial composition in terms of taxa. However, our new analysis shows that gut microbiomes do not necessarily recover all the lost metagenomic sequences. A possible explanation would be that there is a shift in dominant strains or isolates belonging to particular taxa after antibiotic treatment. This metagenomic drift could alter the functional content of microbiomes and affect clinically relevant genetic determinants such as resistance genes.^10,11 In line with the observations of Morowitz and collaborators,¹² our analysis confirms that limiting microbiome studies to taxonomical profiling overlooks the genetic diversity of microbiota.

Figure 2. Comparison of complete assembled metagenomes between A) day 7 vs day 0 and B) day 90 vs day 0 for exposed participants and controls. Horizontal lines indicate the mean of groups. Participants enriched in E. cloacae at day 7 are marked with stars (*).

Figure 3. Comparison of total k-mers content between A) day 7 and day 0 and B) day 90 and day 0. Counts are normalized by the number of reads sequenced per sample. Participants enriched in E. cloacae at day 7 are marked with stars (*).

In our previous work, we created a database of resistance genes (MERGEM database) that can be used for k-mer filtering. ¹³ In this new analysis, we compared metagenomes using MERGEM as a filter for resistance gene k-mers, thereby measuring similarity between gut microbiomes based on sequences potentially associated with antibiotic resistance. The Jaccard similarity index calculated between day 7 and day 0 for resistance gene k-mers in participants with antibiotics (39.0%) was significantly lower than controls (64.9%) (p = 0.023, Welch Two Sample t-test) (Fig. 4A). Difference between exposed (50.1%) and controls (63.4%) for the comparison between day 90 and day 0 is not statistically significant (p = 0.075, Welch Two Sample t-test) (Fig. 4B). At day 7, the number of k-mers previously observed in resistance genes is increased in 8 out of 18 samples (Fig. 5). Notably, participants that had blooms of Enterobacter cloacae complex bacteria had a strong increase in resistance gene-related k-mers compared to controls. The opportunistic pathogen E. cloacae is known as a reservoir of antibiotic resistance genes.^14,15 Our results suggest that a pervasive effect of antibiotics on the gut microbiota could create ecological niches fostering the emergence of multi-resistant opportunistic pathogens, especially in individuals presenting a reduced microbial diversity.

Figure 4. Comparison of resistance gene related k-mers between A) day 7 vs day 0 and B) day 90 vs day 0 for exposed participants and controls. Horizontal lines indicate the mean of groups. Participants enriched in E. cloacae at day 7 are marked with stars (*).

Figure 5. Comparison of k-mers related to antibiotic resistance between A) day 7 and day 0 and B) day 90 and day 0. Counts are normalized by the number of reads sequenced per sample. Participants enriched in E. cloacae at day 7 are marked with stars (*).

As a control experiment, we investigated how the number of reads used to assemble metagenomes affected the Jaccard similarity index for 4 representative participants: 2 controls (C6 and C38) and 2 exposed to antibiotics (E18 and E21). We compared the assemblies generated by subsampling reads from the day 0 samples (with increments of 8 million reads) to the assemblies obtained using all the reads of the day 7 samples. For the 4 samples, we observed that the Jaccard similarity index reached a plateau between 80 and 100 million reads (Fig. 6A). We therefore suggest that a depth of 100 million reads should provide in most cases a consistent approximation of similarity and could limit bias caused by sequencing depth when comparing gut microbiome assemblies using k-mers. Comparing assemblies at a lower sequencing depth using k-mers would generate inconsistent results. Differences in curve shapes between controls and exposed participants (Fig. 6A) might indicate that accurately comparing samples that show little perturbation require more reads than comparing highly affected participants. When looking at how sequencing depth influenced assemblies, we observed that assembly size still continued to grow when more reads were added (Fig. 6B). In contrast, sequencing depth of more than 24 million reads did not seem to affect taxonomical profiling (Fig. 6C).

Figure 6. Influence of sequencing depth on comparison, assembly and profiling of metagenomes for 4 participants (controls 6 and 38, and exposed to antibiotic 18 and 21). A) Comparison of k-mer content of day 0 assemblies with day 7 assemblies using Ray Surveyor. Day 7 assemblies were performed using all the reads sequenced while day 0 assemblies were performed by subsampling reads using 8 million reads increments. Similarity was calculated using the Jaccard similarity index. B) Influence of the number of reads used for assembly on the number of nucleotides assembled into contigs for day 0 samples. Subsampling of reads was in 8 million reads increments. C) Taxonomical profiling at the family rank using Ray Meta 2.3.1 on subsampled reads. Subsampling increments of 8 million nucleotides are indicated by dashes under the plot.

Conclusion

Here, we show that directly comparing the genomic content of microbiomes allows to quantify the extent of microbiome perturbations after antibiotic treatment. This k-mer comparison method takes into account the pan-genome of species and the variability that exists within a taxonomic unit, for example different strains of a species. We can now quantitatively explore how specific strains or genomic elements, such as resistance genes, are modulated by treatment. Since comparisons using k-mers and the Jaccard similarity index are strongly affected by increases in genomic distance,¹⁶ this approach is well suited to analyze time-course experiments or before-and-after states in an individual. Comparison of k-mer content allowed to measure gain or loss of genetic material at the end of the antibiotic treatment and after 3 months of recovery, which could not be observed by taxonomical profiling. However, it remains difficult to answer the critical question of whether a microorganism has disappeared from a population after antibiotic treatment.¹⁷ Indeed, detection of specific taxa is dependent upon sequencing depth, bacterial abundance and microbial diversity.¹⁸

Comparing k-mers to analyze changes in microbiomes following antibiotic stress brings forth the question of how to consider diversity in shotgun microbiome analysis. In most studies, we consider microbial diversity as an index describing population structures that take into account the number of entities measured at a specific taxonomic rank and the evenness of their relative proportions. To extend this concept, researchers are now considering the single nucleotide polymorphisms found in bacterial populations to better understand the effect of environmental conditions or drugs on the strain diversity within metagenomes.^19,20 Unravelling the full richness of microbial genetic diversity in microbiomes will allow a thorough understanding of changes induced by any compounds, drugs or infections.²¹

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Acknowledgment

Computations were performed under the auspices of Calcul Québec and Compute Canada. The operations of Compute Canada are funded by the Canada Foundation for Innovation (CFI), the National Science and Engineering Research Council (NSERC), NanoQuébec and the Fonds Québécois de Recherche sur la Nature et les Technologies (FQRNT).

Funding

This study was financed by the Québec Consortium for Drug Discovery (CQDM) and by the Canada Research Chair Program (JC). FR was supported by a Mitacs post-doctoral fellowship.