Benjamin J. Korry
Damien J. Cabral
Peter Belenky- Department of Molecular Microbiology and Immunology, Division of Biology and Medicine, Brown University, Providence, RI, United States
Antibiotic resistance is a current and expanding threat to the practice of modern medicine. Antibiotic therapy has been shown to perturb the composition of the host microbiome with significant health consequences. In addition, the gut microbiome is known to be a reservoir of antibiotic resistance genes. Work has demonstrated that antibiotics can alter the collection of antibiotic resistance genes within the microbiome through selection and horizontal gene transfer. While antibiotics also have the potential to impact the expression of resistance genes, metagenomic-based pipelines currently lack the ability to detect these shifts. Here, we utilized a dual sequencing approach combining shotgun metagenomics and metatranscriptomics to profile how three antibiotics, amoxicillin, doxycycline, and ciprofloxacin, impact the murine gut resistome at the DNA and RNA level. We found that each antibiotic induced broad, but untargeted impacts on the gene content of the resistome. In contrast, changes in ARG transcript abundance were more targeted to the antibiotic treatment. Doxycycline and amoxicillin induced the expression of tetracycline and beta-lactamase resistance genes, respectively. Furthermore, the increased beta-lactamase resistance gene transcripts could contribute to an observed bloom of Bacteroides thetaiotaomicron during amoxicillin treatment. Based on these findings, we propose that the utilization of a dual sequencing methodology provides a unique capacity to fully understand the response of the resistome to antibiotic perturbation. In particular, the analysis of transcripts reveals that the expression and utilization of resistance genes is far narrower than their abundance at the genomic level would suggest.
Introduction
Antibiotic resistance has emerged as a major threat to human health. In the United States, millions of people suffer from infections caused by antibiotic resistant bacteria, and tens of thousands die as a result (CDC, 2013). Although antibiotic resistance is recognized to be an ancient phenomenon predating the therapeutic use of antibiotics (D’Costa et al., 2011; Bhullar et al., 2012; Santiago-Rodriguez et al., 2015), recent misuse and overuse of antibiotics has led to an increase in the selection for antibiotic resistance genes (ARGs) and has contributed to the spread of infections caused by antibiotic resistant bacteria. Thus, it is important to understand how antibiotic exposure impacts the abundance and expression of resistance genes in the host. Culture-independent methods of profiling entire microbial communities for ARGs have greatly expanded our ability to detect and track resistance elements utilizing high-throughput sequencing techniques (Sommer et al., 2009; Lu et al., 2014; Crofts et al., 2017; Arango-Argoty et al., 2018; Yin et al., 2018). This important development in detection has led to the discovery of ARGs in gut colonizing microbes.
The gut microbiota is now recognized as an important element in human health and disease (Cho and Blaser, 2012; Human Microbiome Project Consortium, 2012; Pflughoeft and Versalovic, 2012; Lynch and Pedersen, 2016), and can serve as a reservoir of antibiotic resistant bacteria (Sommer et al., 2010). Research into the collection of ARGs within the microbiome, termed the “resistome”(D’Costa et al., 2006; Wright, 2007), has begun to explore resistance genes within the gut. Studies have characterized the gut resistome in terms of composition (Sommer et al., 2009, 2010; Hu et al., 2013), life history (Lu et al., 2014; Moore et al., 2015; Parnanen et al., 2018), geographic location (Forslund et al., 2013; Hu et al., 2013; Pehrsson et al., 2016), antibiotic perturbation (Jernberg et al., 2007; Looft et al., 2012), and other factors. In addition, horizontal gene transfer (HGT) of ARGs between bacteria in the gut has been theorized to contribute to the spread of resistance (Shoemaker et al., 2001; Lester et al., 2006; Karami et al., 2007; Smillie et al., 2011; Stecher et al., 2012; Huddleston, 2014). Since the gut microbiome is a reservoir of antibiotic resistance and has the potential to promote HGT of ARGs, it is of particular importance to understand the role of antibiotics in shaping the landscape of antibiotic resistance in this microbial environment.
Antibiotic therapy has been shown to have a dramatic impact on the microbiome, playing a role in gut dysbiosis (Dethlefsen et al., 2008; Dethlefsen and Relman, 2011; Becattini et al., 2016), increasing susceptibility to infection (Ubeda et al., 2010; Buffie et al., 2012), and altering the composition of the gut resistome (Jernberg et al., 2007; Zhang et al., 2013; Palleja et al., 2018). Less is known about how antibiotics promote ARG selection in vivo (Sjolund et al., 2003; Looft et al., 2012), and the impact of antibiotics on the expression of resistance genes in the host microbiota. In response to antibiotic treatment, changes in the resistome may be stochastic, induced by the underlying changes in population structure, or more directed and targeted toward the drug administered. Here we utilize a dual sequencing methodology that employs both metagenomics and metatranscriptomics to reveal broad, untargeted changes in the resistome at the DNA level and a narrower, drug-specific response at the RNA level.
Materials and Methods
Mouse Experiments
In a previous study, we obtained total cecal DNA and RNA from six-week-old female C57BL/6J mice that were treated with amoxicillin (0.1667 mg/mL) for 12 h, or ciprofloxacin (0.0833 mg/mL), or doxycycline hydrochloride (0.067 mg/mL) for 24 h (Cabral et al., 2019). All antibiotic treatments were administered in drinking water, which was provided ad libitum. Based on the estimate that a healthy mouse drinks 150 mL/kg of water per day, these concentrations were selected to administer a dosage of 25 (amoxicillin), 12.5 (ciprofloxacin), or 10 (doxycycline) mg/kg/day (Yang et al., 2017). Control mice were provided with pH-matched water. Each group had four mice that were split into at least two cages per group to account for cage effects. After the 12- or 24-h treatments, mice were sacrificed and cecal contents were collected and stored in DNA/RNA Shield Collection and Lysis Tube from Zymo Research (Irvine, CA, United States). Samples were kept on ice prior to being transferred to −80°C for permanent storage. Mouse experiments were carried out at the Brown University mouse facility with approval from the Institutional Animal Care and Use Committee of Brown University.
Nucleic Acid Extraction and Quantification
DNA and RNA were extracted from cecal samples using the ZymoBIOMICS DNA/RNA Miniprep Kit from Zymo Research (Irvine, CA, United States), eluted in nuclease-free water, and stored at −80°C. Extracted DNA and RNA were quantified using the dsDNA-BR and RNA-HS kits on a QubitTM 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, United States).
Metagenomic and Metatranscriptomic Library Preparation
Metagenomic libraries were prepared using the Ovation® Ultralow System V2 kit from NuGEN (San Carlos, CA, United States). DNA was sheared to a median fragment size of 300 bp using a Covaris S220 High Performance Ultrasonicator (Woburn, MA, United States), and used to prepare metagenomic libraries following the manufacturer’s protocol. Metatranscriptomic libraries were prepared using Ovation® Complete Prokaryotic RNA-seq Library System from NuGEN. Extracted RNA was first treated with DNA rDNase I to remove contaminating DNA. Next, host mRNA and bacterial ribosomal RNA was reduced using the MICROBEnrich and MICROBExpress kits from Invitrogen (Carlsbad, CA, United States). This processed RNA was then used to generate metatranscriptomic libraries following the manufacturer’s protocol with the addition of AnyDeplete probes designed to specifically remove fragments originating from murine osteosarcoma virus, a known source of contamination sequences.
Sequencing
Metagenomic and metatranscriptomic libraries were sequenced on an Illumina HiSeqX using paired-end, 150 bp reads. Sequencing yielded an average of 26,113,145 (± 11,436,616) and 85,599,941 (± 11,674,614) raw reads per metagenomic and metatranscriptomic libraries, respectively. All reads were deposited in the NCBI Short Read Archive under BioProject numbers PRJNA504846 (metagenomics) and PRJNA515074 (metatranscriptomics). This data set was previously published by Cabral et al. (2019). The DNA, RNA, and subsequent libraries and sequencing data were the same as those initially published in Cabral et al. (2019). Here, we reanalyze this data using a different set of pipelines and with a focus on the detection of ARGs.
Processing of Raw Reads
Raw reads from both metagenomic and metatranscriptomic sequencing were processed using the kneadData wrapper script (McIver et al., 2018). Reads were trimmed with Trimmomatic (version 0.36) with SLIDINGWINDOW set at 4:20, MINLEN set at 50, and ILLUMINACLIP: TruSeq3-PE.fa:2:20:10 (Bolger et al., 2014). Sequences from contaminating C57BL/6NJ mouse genome and two murine retroviruses [murine osteosarcoma virus (accession NC_001506.1) and mouse mammary tumor virus (accession NC_001503)] were filtered out using Bowtie2 (Langmead and Salzberg, 2012). In addition to this preprocessing, bacterial ribosomal reads were removed from the datasets using the SILVA 128 database (Pruesse et al., 2007). Based on the principle coordinate analysis (PCoA) of the metatranscriptomic derived resistomes we determined that doxycycline sample 1 was an outlier and it was removed from further analysis (Supplementary Figure S1). Doxycycline sample 1 had roughly 10 times the number of ARG hits relative to all other samples despite it being sequenced to a similar depth (Supplementary Table S1).
Taxonomic Analysis of Metagenomic Reads
Cleaned metagenomic forward reads were classified against a database containing all prokaryotic genomes downloaded from NCBI RefSeq using Kaiju (version 1.7.0) using the MEM run mode and the default cutoffs for E-value and minimum required match length (Menzel et al., 2016) (full relative abundance tables can be found in Supplementary Table S2). The taxonomic output table was analyzed in R (version 3.5.2) using the phyloseq package (version 1.24.2) to calculate alpha and beta diversity metrics (McMurdie and Holmes, 2013). PCoA was performed using the Bray-Curtis Dissimilarity metric (Bray and Curtis, 1957).
Metagenomic Assembly and Binning
The PATRIC web server (3.5.43) (Wattam et al., 2017) was utilized to bin and assign taxonomy to contigs assembled using metaSPAdes within SPAdes (3.13.0) (Nurk et al., 2017).
Antibiotic Resistance Gene Analysis
Processed reads were joined using the fastq-join function of the ea-utils package (Aronesty, 2011). The joined reads were then queried for antibiotic resistance genes using DeepARG (version 1) (Arango-Argoty et al., 2018) using the default settings (0.8 minimum coverage of alignment, E-value cutoff 1e-10, 50% minimum percentage of identity) and excluding “predicted” resistance genes (full counts tables of ARGs and ARG classes can be found in Supplementary Table S3). ARGs at the DNA level from the metagenomic data were identified by running cleaned, merged PE reads through the DeepARG pipeline. On average 45,768 (± 10,213) ARG hits were obtained from an average of 26,281,760 (± 5,410,439) cleaned paired-end (merged forward and reverse read files). Using the same method to identify ARGs at the RNA level from metatranscriptomic sequencing, we found an average of 14,576 (± 6,147) ARG hits from an average of 49,338,997 (± 10,392,548) cleaned paired-end (merged forward and reverse read files) reads.
Statistical Analyses and Figure Generation
Differential abundance of ARGs and ARG classes between treatments and controls was determined using DESeq2 (version 1.20.0) (Love et al., 2014) in R (version 3.5.2) using default parameters. Statistical differences between alpha diversities were calculated in GraphPad Prism (version 8.0) (GraphPad Software, La Jolla, CA, United States). All figures were generated using Prism 8.0, except Supplementary Figure S4 which was generated using the Clustal Omega tool (Sievers and Higgins, 2018) on the EMBL-EBI web server (Madeira et al., 2019).
Results and Discussion
Microbial Diversity
Antibiotic treatment was administered for either 12 h for amoxicillin (beta-lactam) or 24 h for ciprofloxacin (fluoroquinolone) and doxycycline (tetracycline) with untreated time-matched controls. The microbiome of the amoxicillin treated mice displayed a marked reduction in bacterial alpha diversity (p < 0.05, Mann-Whitney U test), however, no change in diversity was observed in mice treated with ciprofloxacin or doxycycline (Figure 1A). Shifts in beta diversity were observed between each treatment and their respective controls indicating that these antibiotics elicit unique changes to the murine gut microbiome at a taxonomic level (p < 0.05, PERMANOVA) (Figure 1B). Perhaps the most drastic taxonomic change was the expansion of Bacteroides thetaiotaomicron in the microbiota of amoxicillin treated mice (Figure 1C and Supplementary Figure S2). Overall, we found that out of the three antibiotics tested, amoxicillin had the most profound impact on the murine cecal microbiome community in terms of diversity and species relative abundance, while ciprofloxacin and doxycycline exhibit less drastic changes. A more detailed description of the taxonomic shifts is detailed in Cabral et al. (2019), while in this study we focus on the ARGs.
Figure 1. Antibiotic treatment has variable impacts on the diversity and taxonomic structure of the microbiome. (A) Shannon diversity index displayed as mean ± SEM (*p < 0.05 Mann-Whitney U test, n = 4). (B) PCoA based on Bray-Curtis beta diversity. (C) Relative abundance of bacterial species displayed as mean ± SEM (n = 4, top 250 most abundant species colored, full relative abundance table available in Supplementary Table S2).
ARG Abundance
While metagenomic analysis is commonly used to characterize the resistome, it can only report the genes found in the community but does not provide information about the actual expression of those genes. In the metagenomic data, we determined the average number of ARGs relative to total reads to calculate the relative abundance of ARGs in the community. There were 1.74E-03 (± 1.62E-04) ARGs per read detected in the metagenomic data (Supplementary Table S1), with an average of 336 (± 14) unique ARGs found in each metagenomic sample (Supplementary Figure S3A). In contrast to the metagenomic data, metatranscriptomics cannot describe the structure of the community, but it can identify the portion of the total genes actively transcribed by the microbiome. In the metatranscriptomic data, there were 3.1E-04 (± 1.67E-04) ARGs per read (Supplementary Table S1), with an average of 161 (± 40) unique ARGs found in each metatranscriptomic sample (Supplementary Figure S3B). This data shows that there are more unique ARGs found in the DNA than in the RNA despite the higher sequencing depth used for metatranscriptomics. This discrepancy could indicate that many of the ARGs encoded in the microbiome are not actively transcribed with or without drug pressure.
Resistome Diversity
Two of the antibiotics examined, amoxicillin and ciprofloxacin, had unique impacts on the taxonomic composition of the microbiome resulting in corresponding shifts in the resistome diversity at the DNA level. Metagenomic data showed that compared to time-matched controls, there was a significant increase in the Shannon diversity of the resistomes with amoxicillin (p < 0.05, Mann-Whitney U test), and a decrease in the Shannon diversity with ciprofloxacin (p < 0.05, Mann-Whitney U test), while doxycycline treatment had no impact (Figure 2A). Analysis of the Bray-Curtis beta diversity revealed significant differences in amoxicillin and ciprofloxacin groups compared to their time-matched controls (p < 0.05, PERMANOVA), but not in the doxycycline treated group (p = 0.2, PERMANOVA) (Figure 2C). However, these shifts in ARG diversity profiles did not necessarily reflect a drug specific selection but rather resulted from the overall shift in microbiome composition. We found that while treatment induced changes in alpha diversity of the resistome at the DNA level, it did not impact resistome alpha diversity at the RNA level (Figure 2B). The lower and more stable alpha diversity of the RNA reads compared to the DNA reads likely stems from the fact that many of the genes detected in the metagenomics are not actively transcribed under vehicle or antibiotic treatment. We also found that amoxicillin and ciprofloxacin did not significantly impact Bray-Curtis ARG diversity at the RNA level (p = 0.128, p = 0.397, PERMANOVA). Additionally, in the metatranscriptomic data, doxycycline did induce a significant Bray-Curtis shift from the 24-h controls (p < 0.05, PERMANOVA) (Figure 2D). This shift in beta-diversity may be driven by the induction of drug targeted ARG transcripts observed in the doxycycline treated samples.
Figure 2. Antibiotics have variable impacts on the diversity and structure of the resistome. (A,B) Shannon diversity index based on resistance gene counts displayed as mean ± SEM for both metagenomic and metatranscriptomic data (*p < 0.05 Mann-Whitney U test). (C,D) PCoA based on Bray-Curtis of resistance gene counts for both metagenomic and metatranscriptomic data.
ARG Class Level Changes in Response to Antibiotics
Antibiotic-induced shifts in ARG classes were determined using DESeq2 and considered significant with a log2FC ≥ 1.5 or ≤ −1.5 and an adjusted p-value < 0.05 (Figure 3A). At the DNA level, we did not find any changes in ARG classes that directly corresponded to antibiotic treatment. For example, the beta-lactam resistance class was not increased with amoxicillin. Instead, we report that the only significant changes in ARG classes were an increase in kasugamycin class ARGs in response to amoxicillin treatment, decreases in the fosmidomycin and trimethoprim classes in response to ciprofloxacin treatment, and a decrease in the fosmidomycin class in response to doxycycline treatment (Figure 3A). As a whole, none of the treatments led to an induction of ARG classes targeted toward the drug administered.
Figure 3. Differential abundance of antibiotic resistance gene classes. Changes in ARG class abundances after antibiotic treatment observed in (A) metagenomic and (B) metatranscriptomic data. The color scale represents log2 fold change (***log2FC ≥ 1.5/ ≤ –1.5 and padj < 0.05) (Full ARG class counts available in Supplementary Table S3).
In contrast to the lack of drug targeted changes at the DNA level, metatranscriptomic sequencing showed an induction of ARG classes targeted to two of the treatments at the RNA level (Figure 3B). Overall, there were a number of significant changes in ARG classes against beta-lactams, fosmidomycin, polymyxin, and trimethoprim in response to amoxicillin treatment, triclosan in response to ciprofloxacin treatment, and fosfomycin, rifampin, and tetracycline in response to doxycycline treatment. Most interestingly, amoxicillin significantly increased the abundance of the beta-lactam ARG class with log2FC 2.67 (padj = 2.50E-04), and doxycycline increased the abundance of the tetracycline ARG class with log2FC 1.81 (padj 3.10E-21) in the RNA (Figure 3B). Thus, in contrast to the metagenomic data, metatranscriptomics shows a significant increase in ARG classes that are targeted to the antibiotic treatment and have the potential to provide a fitness advantage to members of the gut microbiota.
ARG Level Changes in Response to Antibiotics
Results from the differential abundance analysis show that the antibiotics tested have variable impacts on the abundance of AR genes and transcripts. At the metagenomic level, we found a set of differentially abundant genes that appeared general and unrelated to the antibiotic utilized. In contrast, the transcriptional response was much narrower and in the case of amoxicillin and doxycycline, it appears that antibiotic therapy promoted genes directly targeted to the drug utilized. This dichotomy between DNA and RNA level responses could not have been detected without using a dual sequencing approach. Overall, there were fewer differentially abundant ARG transcripts (21 transcripts) found in the metatranscriptomic analysis compared to the number of differentially abundant ARGs (116 genes) in the metagenomic data (Figures 4A–F). This is a reflection of the fewer ARG reads found in the metatranscriptomic data, as well as the more specific response of the microbiome at a transcriptional level compared to the broad metagenomic changes. This is best exemplified by changes in ARGs targeted to antibiotics, specifically amoxicillin and doxycycline.
Figure 4. Differential abundance of antibiotic resistance genes. Changes in ARG abundances after antibiotic treatment observed in metagenomic (A,C,E) and metatranscriptomic data (B,D,F). Bars represent change in gene/transcript abundance after exposure to antibiotics displayed as log2 fold change ± standard error (log2FC ≥ 1.5/ ≤ –1.5 and padj < 0.05) (Full ARG class counts available in Supplementary Table S3).
While the observations made at the ARG class level were fairly broad, the gene level data provided more insights into the direct impact of antibiotic administration on specific ARGs. The differential expression tool DESeq2 was used to analyze the antibiotic-induced changes in AR gene and transcript abundance (Figure 4A). We found 56 significantly elevated or reduced ARGs after amoxicillin treatment. Of these 56, two beta-lactamase genes of interest, cepA and bl2e_cepA, had significant increases in gene abundances of log2FC 4.96 (padj = 4.10E-20) and log2FC 4.73 (padj = 6.72E-16), respectively. In addition to drug targeted genes, we also found increases in a much larger set of untargeted genes (42 genes). It is possible that these changes are the result of taxonomic shifts in bacteria that encode these genes, rather than a direct selection promoted by the induced resistance genes. The beta-lactamase genes increased in the metagenomic data are also increased at the RNA level, highlighted by significantly higher transcript abundances of log2FC 5.12 (padj = 2.40E-05), 4.93, (padj = 2.01E-3), for cepA and bl2e_cepA, respectively (Figure 4B). This may suggest that in response to amoxicillin the community increased transcription of beta-lactamase genes leading to increased bacterial fitness. Conversely, it is also possible that this change merely reflects a bloom in the bacterium encoding these transcripts rather than a direct transcriptional response.
To identify the bacterial origin of the beta-lactamase genes found in our dataset, bacterial genomes were assembled from metagenomic data. Within the B. thetaiotaomicron metagenomically assembled genome (MAG), we identified a region corresponding to a class A beta-lactamase gene with 100% protein sequence homology to a subclass A2 beta-lactamase. Due to its high degree of sequence similarity, this gene likely corresponds to the reads assigned to the cepA and bl2e_cepA genes (Supplementary Figure S4). The relative bloom of B. thetaiotaomicron after amoxicillin treatment may account for the increase in both the cepA and bl2e_cepA gene abundance and transcript level abundance. It is possible that the survival of this taxa during amoxicillin treatment may be promoted by these genes, although we cannot make a definitive conclusion without more evidence. Various laboratory strains and patient isolates of Bacteroidales have been shown to exhibit high levels of resistance to beta-lactams including amoxicillin (Rogers et al., 1994; Nakano et al., 2011; Stentz et al., 2015). Previous research into B. thetaiotaomicron found that this bacterium produces outer membrane vesicles (OMVs) containing cephalosporinase enzymes that protect neighboring bacteria from beta-lactam antibiotics (Stentz et al., 2015). Because B. thetaiotaomicron is a common human commensal (Curtis et al., 2014; Banerjee et al., 2018), this work has interesting implications for the complex microbial environment of the gut microbiome where B. thetaiotaomicron could have a role in modulating antibiotic activity across many taxa through OMV-secreted beta-lactamase enzymes.
Metagenomic data showed that ciprofloxacin treatment induced significant changes in ARG abundance, most notably an increase in the relative abundance of several chloramphenicol, aminoglycoside, and MLS class genes, log2FC > 1.5 (padj < 0.05), and a decrease in the abundance of genes related to multidrug and fosmidomycin resistance log2FC < −1.5 (padj < 0.05) (Figure 4C). No significant increases in fluoroquinolone resistance genes were found in this dataset, however, it should be noted that point mutations conferring resistances were not included in the ARG database used. Thus, fluoroquinolone resistance mutations in gyrA, gyrB, or parC could not be identified from this analysis and might be present in the resistome. No ARG transcripts corresponding to fluoroquinolone resistance were significantly elevated in the ciprofloxacin treated samples. The transcripts of tetC, muxC, and mdtC, all genes encoding efflux system components, were significantly increased by ciprofloxacin treatment (log2FC ≥ 1.5, padj < 0.05) (Figure 4D). Although none of these genes have been shown to directly efflux ciprofloxacin, the fact that fluoroquinolone treatment exclusively increased transcription of efflux type ARGs remains interesting.
At the DNA level, no genes were increased in abundance during doxycycline treatment; however, several were decreased in abundance. Although genes known to confer doxycycline resistance were detected in the metagenomic data (Supplementary Table S3), they were not significantly changed due to doxycycline treatment (Figure 4E). While doxycycline did not increase the abundance of any tetracycline class ARGs in the metagenomic data set, we did detect changes in the metatranscriptomic data. At the RNA level, doxycycline appears to have a targeted effect on the transcript abundance of tetracycline resistance genes with significant increases in the abundance of tet32, tet44, and tetW with log2FCs of 5.06, 4.63, and 4.07, respectively (padj = 1.6E-2, 1.2E-2, 2.2E-2, respectively) (Figure 4E). This distinct difference in gene versus transcript abundance of these tetracycline class ARGs suggests that while short-term treatment with doxycycline may not select for bacteria encoding these resistance genes, it may induce their expression. All three tetracycline resistance genes with increased transcript abundance, tet32, tet44, and tetW have been shown to offer protection to tetracycline antibiotics through ribosomal protection mechanisms (Melville et al., 2001; Stanton et al., 2004; Abril et al., 2010). The increased transcript abundance of several tet genes in response to doxycycline, combined with unchanged levels of these same tet genes in the metagenomic dataset, suggests that their elevated transcriptional activity may be providing protection and enabling the population of tet-carrying bacteria to remain stable during treatment. The doxycycline-induced expression of several tetracycline resistance genes highlights the need for increased transcriptional profiling of ARGs. Relying solely on metagenomics without utilizing metatranscriptomic sequencing, we would have missed this ARG activity that may contribute to bacterial survival during antibiotic pressure.
Conclusion
In this study, we show that short-term antibiotic pressure leads to the differential expression of specific ARGs within the microbiome, but alters the metagenomic landscape in a less targeted way. We use both shotgun metagenomics and metatranscriptomics to profile how three antibiotics, amoxicillin (beta-lactam), doxycycline (tetracycline), and ciprofloxacin (fluoroquinolone), impact the diversity, composition, and transcriptional response of the murine gut resistome. We found that combining these two sequencing methods provides unique perspectives on ARGs in the microbiome that would have been missed by using metagenomics exclusively. For example, we found that at the RNA level a majority of the induced ARGs were targeted against the administered antibiotic, while at the DNA level we found more changes overall, but we did not find a drug-specific pattern. Our results show that both bactericidal and bacteriostatic antibiotic treatment alters the resistome at both the RNA and DNA level, but that these changes may be more specific to the drug administered at the transcript abundance level.
This work highlights the impacts of three antibiotics on the murine cecal resistome as well as the importance of both metagenomic and metatranscriptomic profiling of ARGs. However, due to the current methodology, there are several limitations to this work that must be considered. Experiments were done in mice in a closed mouse facility with reduced exposure to environmental bacteria and the antibiotic exposure period was fairly short. Both of these factors will reduce the opportunity for new ARGs to enter the microbiome and for selection to act on existing ARGs. In addition, due to limitations in strain level identification from current metagenomic methodologies, we are unable to predict whether or not there was selection of specific ARG containing strains following antibiotic treatment. Additionally, the computational pipeline utilized to detect ARGs is unable to identify point mutations and their contributions to the resistome. These aspects limit our ability to detect selective events and the evolution of resistance. Finally, these experiments were conducted with a limited sample size (n = 4), and no samples were collected at time 0. These two factors may limit our ability to detect smaller shifts in ARG levels with sufficient significance and may hinder the detection of baseline changes to the microbiota of the control mice over the course of the experiment.
Despite these limitations, we believe that the dual sequencing approach has real benefits over purely DNA-based approaches. The widespread use of common antibiotics such as those tested in this study contributes to the dissemination of resistance genes in the human population and our microbiomes. However, as demonstrated here, the presence of an ARG in the metagenome may not necessarily indicate that it will be transcribed at baseline or in response to antibiotics. Thus, as we continue to develop strategies to monitor resistance in patient populations it will be important to track both gene presence and expression.
Data Availability Statement
Raw metagenomic and metatranscriptomic reads were deposited in the NCBI Short Read Archive under the BioProject number PRJNA504846.
Ethics Statement
The animal study was reviewed and approved by Institutional Animal Care and Use Committee of Brown University.
Author Contributions
BK, DC, PB contributed to the conception and design of the study. DC performed the experiments and sample processing. BK performed the data analysis and drafted the manuscript. BK, DC, and PB contributed to editing drafts of the manuscript.
Funding
This work was supported by the National Science Foundation through the Graduate Research Fellowship Program under award number 1644760 for BK and DC, by the National Institutes of Health through the National Center for Complementary and Integrative Health award number 1R21A T010366, and by the National Institutes of Health under institutional development awards P20GM121344 and P20GM109035 from the National Institute of General Medical Sciences, which fund the Center for Antimicrobial Resistance and Therapeutic Discovery and the COBRE Center for Computational Biology of Human Disease, respectively. Opinions, interpretations, conclusions, and recommendations are those of the authors and are not necessarily endorsed by the National Science Foundation or the National Institutes of Health.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.00322/full#supplementary-material
References
Abril, C., Brodard, I., and Perreten, V. (2010). Two novel antibiotic resistance genes, tet(44) and ant(6)-Ib, are located within a transferable pathogenicity island in Campylobacter fetus subsp fetus. Antimicrob. Agents Chemother. 54, 3052–3055. doi: 10.1128/aac.00304-10
Arango-Argoty, G., Garner, E., Pruden, A., Heath, L. S., Vikesland, P., and Zhang, L. (2018). DeepARG: a deep learning approach for predicting antibiotic resistance genes from metagenomic data. Microbiome 6:23.
Aronesty, E. (2011). ea-Utils: Command-Line Tools for Processing Biological Sequencing Data. Available at: https://github.com/ExpressionAnalysis/ea-utils (accessed November 16, 2018).
Banerjee, S., Schlaeppi, K., and van der Heijden, M. G. A. (2018). Keystone taxa as drivers of microbiome structure and functioning. Nat. Rev. Microbiol. 16, 567–576. doi: 10.1038/s41579-018-0024-1
Becattini, S., Taur, Y., and Pamer, E. G. (2016). Antibiotic-induced changes in the intestinal microbiota and disease. Trends Mol. Med. 22, 458–478. doi: 10.1016/j.molmed.2016.04.003
Bhullar, K., Waglechner, N., Pawlowski, A., Koteva, K., Banks, E. D., Johnston, M. D., et al. (2012). Antibiotic resistance is prevalent in an isolated cave microbiome. PLoS One 7:e34953. doi: 10.1371/journal.pone.0034953
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Bray, J. R., and Curtis, J. T. (1957). An ordination of the upland forest communities of Southern Wisconsin. Ecol. Monogr. 27, 326–349.
Buffie, C. G., Jarchum, I., Equinda, M., Lipuma, L., Gobourne, A., Viale, A., et al. (2012). Profound alterations of intestinal microbiota following a single dose of clindamycin results in sustained susceptibility to Clostridium difficile-induced colitis. Infect. Immun. 80, 62–73. doi: 10.1128/iai.05496-11
Cabral, D. J., Penumutchu, S., Reinhart, E. M., Zhang, C., Korry, B. J., Wurster, J. I., et al. (2019). Microbial metabolism modulates antibiotic susceptibility within the murine gut microbiome. Cell Metab. 30, 800.e7–823.e7.
CDC (2013). Antibiotic Resistance Threats in the United States. Atlanta, GA: U.S. Department of Health and Human Services.
Cho, I., and Blaser, M. J. (2012). The human microbiome: at the interface of health and disease. Nat. Rev. Genet. 13, 260–270. doi: 10.1038/nrg3182
Crofts, T. S., Gasparrini, A. J., and Dantas, G. (2017). Next-generation approaches to understand and combat the antibiotic resistome. Nat. Rev. Microbiol. 15, 422–434. doi: 10.1038/nrmicro.2017.28
Curtis, M. M., Hu, Z. P., Klimko, C., Narayanan, S., Deberardinis, R., and Sperandio, V. (2014). The gut commensal bacteroides thetaiotaomicron exacerbates enteric infection through modification of the metabolic landscape. Cell Host Microbe 16, 759–769. doi: 10.1016/j.chom.2014.11.005
D’Costa, V. M., King, C. E., Kalan, L., Morar, M., Sung, W. W., Schwarz, C., et al. (2011). Antibiotic resistance is ancient. Nature 477, 457–461.
D’Costa, V. M., McGrann, K. M., Hughes, D. W., and Wright, G. D. (2006). Sampling the antibiotic resistome. Science 311, 374–377. doi: 10.1126/science.1120800
Dethlefsen, L., Huse, S., Sogin, M. L., and Relman, D. A. (2008). The pervasive effects of an antibiotic on the human gut microbiota, as revealed by deep 16S rRNA sequencing. PLoS Biol. 6:e280. doi: 10.1371/journal.pbio.0060280
Dethlefsen, L., and Relman, D. A. (2011). Incomplete recovery and individualized responses of the human distal gut microbiota to repeated antibiotic perturbation. Proc. Natl. Acad. Sci. U.S.A. 108(Suppl. 1), 4554–4561. doi: 10.1073/pnas.1000087107
Forslund, K., Sunagawa, S., Kultima, J. R., Mende, D. R., Arumugam, M., Typas, A., et al. (2013). Country-specific antibiotic use practices impact the human gut resistome. Genome Res. 23, 1163–1169. doi: 10.1101/gr.155465.113
Hu, Y., Yang, X., Qin, J., Lu, N., Cheng, G., Wu, N., et al. (2013). Metagenome-wide analysis of antibiotic resistance genes in a large cohort of human gut microbiota. Nat. Commun. 4:2151.
Huddleston, J. R. (2014). Horizontal gene transfer in the human gastrointestinal tract: potential spread of antibiotic resistance genes. Infect. Drug Resist. 7, 167–176.
Human Microbiome Project Consortium (2012). Structure, function and diversity of the healthy human microbiome. Nature 486, 207–214. doi: 10.1038/nature11234
Jernberg, C., Lofmark, S., Edlund, C., and Jansson, J. K. (2007). Long-term ecological impacts of antibiotic administration on the human intestinal microbiota. ISME J. 1, 56–66. doi: 10.1038/ismej.2007.3
Karami, N., Martner, A., Enne, V. I., Swerkersson, S., Adlerberth, I., and Wold, A. E. (2007). Transfer of an ampicillin resistance gene between two Escherichia coli strains in the bowel microbiota of an infant treated with antibiotics. J. Antimicrob. Chemother. 60, 1142–1145. doi: 10.1093/jac/dkm327
Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9:357. doi: 10.1038/nmeth.1923
Lester, C. H., Frimodt-Moller, N., Sorensen, T. L., Monnet, D. L., and Hammerum, A. M. (2006). In vivo transfer of the vanA resistance gene from an Enterococcus faecium isolate of animal origin to an E. faecium isolate of human origin in the intestines of human volunteers. Antimicrob. Agents Chemother. 50, 596–599. doi: 10.1128/aac.50.2.596-599.2006
Looft, T., Johnson, T. A., Allen, H. K., Bayles, D. O., Alt, D. P., Stedtfeld, R. D., et al. (2012). In-feed antibiotic effects on the swine intestinal microbiome. Proc. Natl. Acad. Sci. U.S.A. 109, 1691–1696. doi: 10.1073/pnas.1120238109
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550.
Lu, N., Hu, Y., Zhu, L., Yang, X., Yin, Y., Lei, F., et al. (2014). DNA microarray analysis reveals that antibiotic resistance-gene diversity in human gut microbiota is age related. Sci. Rep. 4:4302.
Lynch, S. V., and Pedersen, O. (2016). The human intestinal microbiome in health and disease. N. Engl. J. Med. 375, 2369–2379. doi: 10.1056/nejmra1600266
Madeira, F., Park, Y. M., Lee, J., Buso, N., Gur, T., Madhusoodanan, N., et al. (2019). The EMBL-EBI search and sequence analysis tools APIs in 2019. Nucleic Acids Res. 47, W636–W641.
McIver, L. J., Abu-Ali, G., Franzosa, E. A., Schwager, R., Morgan, X. C., Waldron, L., et al. (2018). bioBakery: a meta’omic analysis environment. Bioinformatics 34, 1235–1237. doi: 10.1093/bioinformatics/btx754
McMurdie, J., and Holmes, S. (2013). phyloseq: an r package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8:e61217. doi: 10.1371/journal.pone.0061217
Melville, C. M., Scott, K. P., Mercer, D. K., and Flint, H. J. (2001). Novel tetracycline resistance gene, tet(32), in the Clostridium-related human colonic anaerobe K10 and its transmission in vitro to the rumen anaerobe Butyrivibrio fibrisolvens. Antimicrob. Agents Chemother. 45, 3246–3249. doi: 10.1128/aac.45.11.3246-3249.2001
Menzel, P., Ng, K. L., and Krogh, A. (2016). Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat. Commun. 7:11257.
Moore, A. M., Ahmadi, S., Patel, S., Gibson, M. K., Wang, B., Ndao, M. I., et al. (2015). Gut resistome development in healthy twin pairs in the first year of life. Microbiome 3:27.
Nakano, V., Nascimento E Silva, A., Merino, V. R., Wexler, H. M., and Avila-Campos, M. J. (2011). Antimicrobial resistance and prevalence of resistance genes in intestinal Bacteroidales strains. Clinics 66, 543–547. doi: 10.1590/s1807-59322011000400004
Nurk, S., Meleshko, D., Korobeynikov, A., and Pevzner, A. (2017). metaSPAdes: a new versatile metagenomic assembler. Genome Res. 27, 824–834. doi: 10.1101/gr.213959.116
Palleja, A., Mikkelsen, K. H., Forslund, S. K., Kashani, A., Allin, K. H., Nielsen, T., et al. (2018). Recovery of gut microbiota of healthy adults following antibiotic exposure. Nat. Microbiol. 3, 1255–1265. doi: 10.1038/s41564-018-0257-9
Parnanen, K., Karkman, A., Hultman, J., Lyra, C., Bengtsson-Palme, J., Larsson, D. G. J., et al. (2018). Maternal gut and breast milk microbiota affect infant gut antibiotic resistome and mobile genetic elements. Nat. Commun. 9:3891.
Pehrsson, E. C., Tsukayama, P., Patel, S., Mejia-Bautista, M., Sosa-Soto, G., Navarrete, K. M., et al. (2016). Interconnected microbiomes and resistomes in low-income human habitats. Nature 533, 212–216. doi: 10.1038/nature17672
Pflughoeft, K. J., and Versalovic, J. (2012). Human microbiome in health and disease. Annu. Rev. Pathol. 7, 99–122.
Pruesse, E., Quast, C., Knittel, K., Fuchs, B. M., Ludwig, W. G., Peplies, J., et al. (2007). SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 35, 7188–7196. doi: 10.1093/nar/gkm864
Rogers, M. B., Bennett, T. K., Payne, C. M., and Smith, C. J. (1994). Insertional activation of cepA leads to high-level beta-lactamase expression in Bacteroides fragilis clinical isolates. J. Bacteriol. 176, 4376–4384. doi: 10.1128/jb.176.14.4376-4384.1994
Santiago-Rodriguez, T. M., Fornaciari, G., Luciani, S., Dowd, S. E., Toranzos, G. A., Marota, I., et al. (2015). Gut microbiome of an 11th Century A.D. pre-Columbian andean mummy. PLoS One 10:e0138135. doi: 10.1371/journal.pone.0138135
Shoemaker, N. B., Vlamakis, H., Hayes, K., and Salyers, A. A. (2001). Evidence for extensive resistance gene transfer among Bacteroides spand among Bacteroides and other genera in the human colon. Appl. Environ. Microbiol. 67, 561–568. doi: 10.1128/aem.67.2.561-568.2001
Sievers, F., and Higgins, D. G. (2018). Clustal Omega for making accurate alignments of many protein sequences. Protein Sci. 27, 135–145. doi: 10.1002/pro.3290
Sjolund, M., Wreiber, K., Andersson, D. I., Blaser, M. J., and Engstrand, L. (2003). Long-term persistence of resistant Enterococcus species after antibiotics to eradicate Helicobacter pylori. Ann. Intern. Med. 139, 483–487.
Smillie, C. S., Smith, M. B., Friedman, J., Cordero, O. X., David, L. A., and Alm, E. J. (2011). Ecology drives a global network of gene exchange connecting the human microbiome. Nature 480, 241–244. doi: 10.1038/nature10571
Sommer, M. O., Church, G. M., and Dantas, G. (2010). The human microbiome harbors a diverse reservoir of antibiotic resistance genes. Virulence 1, 299–303. doi: 10.4161/viru.1.4.12010
Sommer, M. O. A., Dantas, G., and Church, G. M. (2009). Functional characterization of the antibiotic resistance reservoir in the human microflora. Science 325, 1128–1131. doi: 10.1126/science.1176950
Stanton, T. B., McDowall, J. S., and Rasmussen, M. A. (2004). Diverse tetracycline resistance genotypes of Megasphaera elsdenii strains selectively cultured from swine feces. Appl. Environ. Microbiol. 70, 3754–3757. doi: 10.1128/aem.70.6.3754-3757.2004
Stecher, B., Denzler, R., Maier, L., Bernet, F., Sanders, M. J., Pickard, D. J., et al. (2012). Gut inflammation can boost horizontal gene transfer between pathogenic and commensal Enterobacteriaceae. Proc. Natl. Acad. Sci. U.S.A. 109, 1269–1274. doi: 10.1073/pnas.1113246109
Stentz, R., Horn, N., Cross, K., Salt, L., Brearley, C., Livermore, D. M., et al. (2015). Cephalosporinases associated with outer membrane vesicles released by Bacteroides spprotect gut pathogens and commensals against beta-lactam antibiotics. J. Antimicrob. Chemother. 70, 701–709. doi: 10.1093/jac/dku466
Ubeda, C., Taur, Y., Jenq, R. R., Equinda, M. J., Son, T., Samstein, M., et al. (2010). Vancomycin-resistant Enterococcus domination of intestinal microbiota is enabled by antibiotic treatment in mice and precedes bloodstream invasion in humans. J. Clin. Invest. 120, 4332–4341. doi: 10.1172/jci43918
Wattam, A. R., Davis, J. J., Assaf, R., Boisvert, S., Brettin, T., Bun, C., et al. (2017). Improvements to PATRIC, the all-bacterial bioinformatics database and analysis resource center. Nucleic Acids Res. 45, D535–D542.
Wright, G. D. (2007). The antibiotic resistome: the nexus of chemical and genetic diversity. Nat. Rev. Microbiol. 5, 175–186. doi: 10.1038/nrmicro1614
Yang, J. H., Bhargava, P., McCloskey, D., Mao, N., Palsson, B. O., and Collins, J. J. (2017). Antibiotic-induced changes to the host metabolic environment inhibit drug efficacy and alter immune function. Cell Host Microbe 22, 757.e3–765.e3.
Yin, X., Jiang, X. T., Chai, B., Li, L., Yang, Y., Cole, J. R., et al. (2018). ARGs-OAP v2.0 with an expanded SARG database and hidden markov models for enhancement characterization and quantification of antibiotic resistance genes in environmental metagenomes. Bioinformatics 34, 2263–2270. doi: 10.1093/bioinformatics/bty053
Keywords: antibiotic resistance genes, microbiome, resistome, antibiotics, metagenomics, metatranscriptomics
Citation: Korry BJ, Cabral DJ and Belenky P (2020) Metatranscriptomics Reveals Antibiotic-Induced Resistance Gene Expression in the Murine Gut Microbiota. Front. Microbiol. 11:322. doi: 10.3389/fmicb.2020.00322
Received: 25 November 2019; Accepted: 13 February 2020;
Published: 06 March 2020.
Edited by:
Jose L. Martinez, Spanish National Research Council, SpainReviewed by:
Felipe Lira, Institut National de la Recherche Agronomique (INRA), FranceSilvia García Cobos, Carlos III Health Institute, Spain
Copyright © 2020 Korry, Cabral and Belenky. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Peter Belenky, cGV0ZXJfYmVsZW5reUBicm93bi5lZHU=