Abstract
The β-lactams are the largest group of clinically applied antibiotics, and resistance to these is primarily associated with β-lactamases. There is increasing understanding that these enzymes are ubiquitous in natural environments and henceforth, elucidating the global diversity, distribution, and mobility of β-lactamase-encoding genes is crucial for holistically understanding resistance to these antibiotics. In this study, we screened 232 shotgun metagenomes from ten different environments against a custom-designed β-lactamase database, and subsequently analyzed β-lactamase homologs with a suite of bioinformatic platforms including cluster and network analyses. Three interrelated β-lactamase clusters encompassed all of the human and bovine feces metagenomes, while β-lactamases from soil, freshwater, glacier, marine, and wastewater grouped within a separate “environmental” cluster that displayed high levels of inter-network connectivity. Interestingly, almost no connectivity occurred between the “feces” and “environmental” clusters. We attributed this in part to the divergence in microbial community composition (dominance of Bacteroidetes and Firmicutes vs. Proteobacteria, respectively). The β-lactamase diversity in the “environmental” cluster was significantly higher than in human and bovine feces microbiomes. Several class A, B, C, and D β-lactamase homologs (blaCTX-M, blaKPC, blaGES) were ubiquitous in the “environmental” cluster, whereas bovine and human feces metagenomes were dominated by class A (primarily cfxA) β-lactamases. Collectively, this study highlights the ubiquitous presence and broad diversity of β-lactamase gene precursors in non-clinical environments. Furthermore, it suggests that horizontal transfer of β-lactamases to human-associated bacteria may be more plausible from animals than from terrestrial and aquatic microbes, seemingly due to phylogenetic similarities.
Introduction
The β-lactam antibiotics are the most widely used clinically applied antibiotics in the world (Garau et al., 2005; European Centre for Disease Prevention and Control [ECDC], 2009; European Surveillance of Antimicrobial Consumption [ESAC], 2009) together with cephalosporins and broad-spectrum penicillins constituting approximately 55% of globally consumed antibiotics (Van Boeckel et al., 2014). Nonetheless, the success of antibiotic therapy is compromised by the development of resistance in many important clinical pathogens (Davies and Davies, 2010). The major cause of resistance to β-lactam antibiotics are β-lactamases, which deactivate antibacterial properties by hydrolysis (Davies, 1994). β-lactamases are found in a broad range of bacteria and include multiple genes that are classified in groups or classes based on functional characteristics and sequence similarities, respectively (Medeiros, 1997; D’Costa et al., 2007). In this context, β-lactamases with similar nucleotide sequences often have similar functional properties, but low sequence similarity is not an indication of different functional properties since genes with low similarity may have similar tertiary structure (Kelly et al., 1982; D’Costa et al., 2007). Furthermore, these genes are often associated with mobile genetic elements (MGEs) that facilitate their transfer within and between bacterial taxa (El Salabi et al., 2013).
Despite evidence that antibiotic resistance (AR) is ubiquitous in non-clinical microbiomes, most studies have assessed the incidence and diversity of β-lactamases in pathogens (Rice, 2001; Paterson and Bonomo, 2005; Ramphal and Ambrose, 2006; Shahid, 2014; Sullivan and Schaus, 2015), and only a few have explored the diversity and distribution of these genes in non-clinical environments. Previous studies that screened β-lactamases in the environment have applied culture-dependent approaches (Villedieu et al., 2003; Gijón et al., 2012; Walsh and Duffy, 2013), culture-independent qPCR-based analyses that target particular antibiotic resistance genes (ARGs) (Pallecchi et al., 2007; Donato et al., 2010; Popowska et al., 2012) and functional metagenomics (Allen et al., 2009; Torres-Cortés et al., 2011; Amos et al., 2014; Forsberg et al., 2014; Su et al., 2014) that provide evidence of β-lactamases activity for novel genes. Recent studies have also applied shotgun metagenomic approaches to assess the scope and abundance of ARGs in natural environments (Yang et al., 2013; Nesme et al., 2014; Li et al., 2015; Fondi et al., 2016). Unlike culture-based methods and functional metagenomics, shotgun metagenomics relies on bioinformatic annotations and therefore, does not provide concrete evidence that β-lactamases are actually active. Nonetheless, this approach provides a global perspective on gene abundance and diversity in targeted microbiomes, at a depth unattainable with the other methods.
Elucidating the global distribution of β-lactamases has great importance. Anthropogenically derived genes can be disseminated to natural ecosystems through sewage, treated wastewater, animal feces, and other sources. Furthermore, β-lactamases have a long evolutionary history in natural environments (Hall and Barlow, 2004; Garau et al., 2005; Aminov, 2009; D’Costa et al., 2011), and there is evidence suggesting that these are reservoirs of novel β-lactamase gene precursors that can be transferred to pathogenic bacteria and contribute to clinically associated AR (D’Costa et al., 2007; Berglund, 2015; Versluis et al., 2015). Despite recent studies that explored resistomes in a range of distinct environments (Bhullar et al., 2012; Cristóbal-Azkarate et al., 2014; Nesme et al., 2014; Czekalski et al., 2015; Fondi et al., 2016), there is a lack of collective knowledge regarding the global diversity and distribution of β-lactamase gene homologs. Although circumstantial evidence implies that β-lactamases have mobilized from natural microbial communities to human-associated microbiomes (Forsberg et al., 2012), the extent of this phenomenon, and the overall dissemination of β-lactamases across ecological and geographic boundaries, is currently unknown.
In this study, we screened 232 metagenomes from ten distinct environments against a collated non-redundant β-lactamase database, and subsequently applied different analyses in order to elucidate the diversity, distribution, and abundance of β-lactamase homologs across environmental and geographic compartments.
Materials and Methods
Metagenomic Data Sets
Illumina shotgun metagenomic sequences used in this study were downloaded from two repositories, MG-RAST1 and EBI METAGENOMICS2, after quality control step performed by each repository. Altogether, we retrieved 232 metagenomes from 27 different projects encompassing 4.7 billion sequences (Supplementary Table S1). The metagenomes were associated with ten environments including soil (non-agricultural and agricultural), freshwater, marine, glaciers, human (from healthy people) and bovine feces, bovine rumen, aerobic and anaerobic compartments of municipal wastewater treatment facilities and anaerobic reactors for industrial food waste digestion.
Construction of a β-Lactamase Database
A consolidated β-lactamase genes database coined EX-B, integrated four clinically important publically available databases: The Lahey β-lactamase database, The Lactamase Engineering Database (LACED), The Comprehensive Antibiotic Resistance Database (CARD) and The Pasteur Institute’s OXY, OKP, and LEN protein variation databases (Gatica et al., 2016). Sequences deposited into EX-B were aligned using BioEdit 7.2.5 software, generating 1566 non-redundant β-lactamase sequences. The EX-B database is available for download3. The short reads generated by Illumina sequencing prevent characterization of specific β-lactamase variants, and therefore we collated variants associated with specific families (i.e., all of the blaOXA homologs), resulting in 66 unique β-lactamase gene groups.
Identification and Characterization of β-Lactamase Gene Homologs in Environmental Metagenomes
Characterization of β-lactamase homologs in the targeted metagenomes was achieved by comparing generated reads to the EX-B database using blastx (Gish and States, 1993), with an identity cutoff greater than 50%, a bit score higher than 30 and an e-value lower than 1e-4. Although relatively low, we chose this cutoff because environmental β-lactamase homologs can be diverse relative to clinically characterized β-lactamase genes that dominate public databases. The selected cutoff was validated by contrast, using blastx, randomly selected reads coming from different environments against NCBI non-redundant database.
To reduce overlapping results we only registered the best hit for each read. Normalization of sample read abundance enabled comparison of samples with different sequencing depth; and variance analyses (ANOVA with Tukey and non-parametric Wilcoxon tests for unequal number of observations) were performed to identify differences in the relative abundance and β-lactamase gene diversity across the environments analyzed.
We specifically targeted four genes (blaOXA, blaTEM, blaCTX-M, and blaGES) to determine the level of identity of the reads from the corresponding environments to the EX-B database. Five specific categories were defined based on sequence similarity (50–59%; 60–69%; 70–79%; 80–89%; and 90–100%).
Ordination of β-Lactamase Genes Assemblages in the Targeted Metagenomes
Identified β-lactamase hits from the targeted metagenomes were analyzed in PC-ORD 5.0 (McCune and Mefford, 2011). Initially, the hits were relativized by weighting by ubiquity and transformed by square root. Outliers were identified and removed, and a Bray-Curtis distance matrix using variance regression as an endpoint selection method and Euclidean residual distances to produce scores ranging from 0 to 1 was constructed for use in downstream analysis.
In addition, indicator species analysis was performed in order to pinpoint representative β-lactamases in the targeted environments. This approach is based on a Dufrene and Legendre method that combines information on species abundance and their faithfulness of occurrence in a particular environment. The analysis produces scores from zero (no indication) to 100 (perfect indication). Subsequently, a randomization test (Monte Carlo) based on 4999 runs evaluates the statistical significance of the maximum indicator value recorded for a given species, indicating if the probability to find a given species in a given group is high or low (McCune and Mefford, 2011).
Filtered and normalized hits were analyzed in PC-ORD 5.0 using non-metric multidimensional scaling (NMDS) with a post hoc multi-response permutation process (MRPP) to test for similarities, and subsequently with an indicator species analysis to pinpoint β-lactamases in the targeted environments. In addition, Bray-Curtis distance matrices were generated to construct β-lactamase gene networks using EDENetwork 2.18 (Kivelä et al., 2014), after removing β-lactamase hits that could not be characterized to class level. Finally, graphical network visualization and statistical network tests were performed using Cytoscape 3.4.0 (Shannon, 2003).
Determination of Phylogenetic Distribution in the Targeted Environments
The bacterial community structure of selected environments (agricultural and non-agricultural soils, marine, freshwater, wastewater, human and bovine feces and bovine rumen) was determined by randomly selecting the 16S rRNA phylum-level distributions from 47 metagenomes (Supplementary Table S2), then the average relative abundance per environment was obtained and compared.
Results
Richness of β-Lactamases Homologs in the Environment
The abundance of β-lactamase gene homologs in the ten targeted environments was on average 0.01% in soil [n = 80 metagenomes, standard deviation (STDEV) = 0.002], 0.0068% in glacier (n = 10, STDEV = 0.001), 0.0046% in freshwater (n = 11, STDEV = 0.001), 0.002% in marine (n = 22, STDEV = 0.003), 0.0121 in human feces (n = 63, STDEV = 0.007), 0.0049% in bovine (n = 27, including feces and rumen samples, overall STDEV = 0.004), and 0.0092% in wastewater treatment plant (n = 19), including municipal wastewater (activated sludge and anaerobic digesters) and food waste anaerobic digesters, overall STDEV = 0.003) metagenomes (Supplementary Figure S1 and Supplementary Table S1). Metagenome reads were designated as specific β-lactamase homologs by comparing them to the EX-B database using the blastx algorithm as described in the materials and methods section. Collectively, soil (including agricultural and non-agricultural soils) and glacier metagenomes showed the highest diversity (37 different β-lactamase gene types on average), followed by wastewater metagenomes (30 and 32 different β-lactamase gene types in municipal wastewater treatment facilities and anaerobic digesters, respectively; Figure 1).
Figure 1
Analysis of the detected β-lactamase homologs at class level revealed that class A were dominant in bovine and human feces samples (more than 70% in each case); class B were profuse in soils, bovine rumen, marine, and wastewater treatment plants (40–60% in each environment); class C were abundant in soils, fresh water and bovine rumen environments (approx. 20% in each environment); and class D were copious in glaciers, wastewater and freshwater environments (approx. 15–30%; Supplementary Figure S2). The relative abundance of the five most dominant β-lactamase gene homologs in each environment was also assessed (Figure 2). Homologs of cfxA were dominant in human and bovine feces, whereas blaOXA, blaLRA and L1 metalo-β-lactamase homologs were dominant in natural environments and wastewater treatment facilities. Overall, there were no differences in distribution of the five dominant β-lactamase homologs in non-agricultural vs. agricultural soils, and between freshwater and soil samples (agricultural and non-agricultural soils) (p > 0.05, Wilcoxon test).
Figure 2
Indicator species analysis was applied to pinpoint “core” β-lactamase homologs in selected environments. Four ecosystems with varying inferred levels of “anthropogenic” impact were targeted (human feces, wastewater, agricultural soils, and less anthropogenically impacted environments, which included non-agricultural soils, glaciers, marine, and freshwater samples). There was an inverse relation between the extent of “anthropogenicity,” and the number of indicator species observed. In total, 25 β-lactamase family homologs showed a high probability of occurrence in the less anthropogenically impacted environments, and only four (blaEBRcfxA, HGI, and mecR) exhibited a high probability of occurrence in the human fecal samples (Table 1).
Table 1
| Environment | High probability of β -lactamase homolog |
|---|---|
| Human feces | blaEBR, cfxA, HGI, mecR |
| Wastewater treatment plants | bla1, blaFEZ, blaGOB, blaNPS, blaSIM, blaTEM, blaTOHO, blaZ, blm, LCR |
| Agricultural soils | AmpC, blaACT, blaB, blaLEN, blaLRA, blaMOX, blaNDM, blaOCH, blaSHV, blaV IM, CepS, THINB |
| Less anthropogenically impacted (include non-agricultural soils, glaciers, marine, and fresh water samples) | blaACC, blaCARB, blaCMY, blaCTX-M, blaDHA, blaFONA, blaGES, blaIND, blaJOHN, blaKHM, blaKLUA, blaKPC, blaMIR, blaOKP, blaOXA, blaOXY, blaPDC, blaPER, blaSME, blaSPM, blaTUS, blaV EB, CphA, L1, SRT |
Faithfulness of β-lactamase homolog occurrence based on indicator species analysis in selected metagenomes with varying levels of anthropogenic impact.
Inter- and Intra-Environmental β-Lactamase Associations
We constructed distance matrices and gene networks from all of the metagenome derived β-lactamase hits, and exclusively from well-defined hits characterized as specific β-lactamase families (Supplementary Figure S3), this means that hits defined only as β-lactamase or class β-lactamase were discarded from analysis. In both cases, the topography of the generated β-lactamase gene networks formed eight clusters, with a high level of connectivity within the individual clusters, but with low connectivity between the clusters (Supplementary Tables S3, S4). Cluster I collated all of the soil, freshwater, marine, glacier, and municipal wastewater metagenome-derived β-lactamase homologs; cluster II harbored the food waste anaerobic digester samples; cluster III contained the bovine feces samples; cluster IV and V contained human feces samples (from Europe and the United States, respectively); cluster VI and VII contained bovine rumen samples; and cluster VIII contained freshwater samples. Color coding of the metagenomes by both environmental (Figure 3A) and geographic (Figure 3B) origin revealed high connectivity between metagenomes from similar environments, and limited or no connectivity between samples defined by geographic origin. In some cases such as in the human feces samples, geographic influence occurred within individual environments. Cluster IV encompassed human feces metagenomes from Northern and Central Europe (n = 44), whereas closely associated cluster V contained all of the metagenomes from the United States (n = 15). This geographic differentiation should, however, be taken with caution considering the potential impact of additional differentiating factors such as age and gender. In order to highlight potential interactions between β-lactamase homologs from the natural ecosystems encompassed in Cluster I, we applied several network indices including closeness centrality, radiality, topological coefficients, neighborhood connectivity, and number of directed edges (Supplementary Tables S5, S6). In general, Cluster I exhibited high values of several measures of network centrality confirming the close relationship of the β-lactamase homologs in these environments. Specifically, most of the non-agricultural soil metagenomes (und_1, 4, 12, 22, 25, 27, 28, 30, 31, 55, and 56) were centrally situated in this cluster and exhibited the highest level of connectivity along with a glacier sample (glac_2), and two freshwater samples (fresh_8 and 17). In contrast, connectivity was generally lower in the agricultural soil samples, and these were less centrally situated within the cluster. Evidence suggested that the observed divergence in the agricultural soil samples was strongly associated with plant species, but this could not be statistically validated.
Figure 3
We specifically assessed the microbiomes (at phylum level based on 16S rRNA gene analysis) of eight selected environments (Supplementary Figure S4) to determine potential correlations between β-lactamase composition and bacterial community structure. Collectively, this phylogenetic characterization revealed two groups. “Anaerobic environments,” encompassing human and animal feces and rumen samples that contained high abundances of Firmicutes (relative abundance of 33.6–50.4%) and Bacteroidetes (19.7–24.7%); and “aerobic environments,” dominated by Proteobacteria (22.7–49.4%). Based on this classification, the diversity of β-lactamases was lower in anaerobic environments (although this was not always statistically significant), but the β-lactamase richness across environments did not show differences between the aerobic and anaerobic environments.
Similarity of Environmental β-Lactamase Homologs to Clinical-Associated β-Lactamases
To elucidate the level of similarity between environmental homologs and β-lactamase genes from pathogenic bacteria, we specifically targeted four prevalent clinically associated genes: blaCTX-M, blaGES, blaOXA, and blaTEM, in the metagenomes of the targeted environments. This was based on the notion that the β-lactamase sequences in the EX-B database used in this analysis are predominantly associated with clinically derived bacteria. Generated BLAST hits were grouped in five categories according to the level of identity of the environmental β-lactamase homologs to the database sequences (50–59%; 60–69%; 70–79%; 80–89%; and 90–100%; Figure 4). For the most part, the natural environments showed relatively low levels of homology to the four β-lactamases. In contrast, β-lactamase homologs from municipal wastewater treatment plants, anaerobic food and wastewater digesters and human feces metagenomes displayed higher levels of similarity, suggesting that these might be more “anthropogenically impacted.” Although blaOXA was the most abundant β-lactamase in all samples, blaTEM showed higher levels of similarity to the database sequences, especially in the marine and human feces metagenomes. Interestingly, almost all of the highly similar (90–100%) blaTEM hits in the human feces were related to one specific metagenome (H_gut 8, 1196 hits), suggesting that this sample may be outside the norm and potentially associated with dysbiosis due to antibiotic therapy. Collectively, most of the hits from all of the targeted environments showed less than 70% homology to these four targeted β-lactamases, although high levels of database identity (between 90 and 100%) were observed in some of the targeted environments such as human feces (blaTEM and blaOXA), oceans (blaTEM), and wastewater (blaGES) metagenomes.
Figure 4
We plotted and visualized the abundance and distribution of the four selected β-lactamase genes across the targeted environments using NMDS ordination (Figure 5). Homologs associated with blaOXA were the most widely distributed, with ubiquitous presence in soil, fresh water, marine, and human feces metagenomes. Conversely, blaCTX-M homologs were less abundant, and primarily identified in soils, while blaGES homologs were profuse in soils and in wastewater environments. Interestingly, despite the copiousness of blaTEM in clinical environments, it was not highly distributed across the analyzed environments.
Figure 5
Finally, to support the application of the 50% similarity cutoff, we applied the blastx algorithm to randomly selected reads (n = 40) coming from different environments analyzed in this study, against the NCBI non-redundant gene database (Supplementary Table S7). The results indicated that in only five instances, the first hit was not a β-lactamase, and of these, only one query was not characterized as a β-lactamase in the first ten hits obtained. In general, the EX-B database was more restrictive than NCBI database in 31 cases.
Discussion
Understanding ARG dynamics on a global scale is highly complex. As stipulated in the resistome hypothesis (D’Costa et al., 2006), the long evolutionary history of ARGs makes soil and other natural ecosystems highly diverse reservoirs of ARG precursors that can potentially be transferred to human pathogens (Hall and Barlow, 2004; Garau et al., 2005; Baltz, 2008; Aminov, 2009). Concurrently, dissemination of ARGs from anthropogenic sources to natural environments (Czekalski et al., 2015) can affect the structure of these natural resistomes and contribute to the global scope of AR. In order to enhance understanding of the environmental dimensions of AR and elucidate relationships between natural and anthropogenically associated resistomes, we evaluated the scope, diversity and distribution of β-lactamase gene homologs on a global scale by targeting a multitude of shotgun metagenomes from ten distinct environments. We applied relatively modest cutoff criteria (identity ≥ 50%; bit score ≥ 30; e-value ≤ 1e-4) for defining specific β-lactamases. This is lower than some previously conducted studies (i.e., Yang et al., 2016 used identity ≥ 60%; e-value ≤ 1e-7), but similar to other studies that used lower cutoffs to pinpoint ARGs in metagenomics reads (i.e., Nesme et al., 2014 used identity ≥ 35%; bit score ≥ 60). Application of lower stringency may capture many β-lactamase gene homologs such as esterases (Wagner et al., 2002), proteases (Delfosse et al., 2009), and transpeptidases (Massova and Mobashery, 1998); and therefore, our analyses cannot determine that the pinpointed homologs are actually functional β-lactamases. Conversely, the relation between environmental ARG precursors and homologous clinical genes that dominate public databases is often low (Forsberg et al., 2014) and therefore, by applying highly stringent cutoffs, we may omit a substantial fraction of the environmental resistome. The selected cutoff was validated by contrast, using blastx, randomly selected reads against NCBI non-redundant database; thus, the results indicated that even when a low similarity was observed (50%), the reads were more close to ARGs than to any other type of gene. Even more, usually, reads contrasted against NCBI database, presented a higher percent of identity to β-lactamases than the observed using EX-B database.
The β-lactamase homologs pinpointed in this study distinctly grouped into four mega-clusters (human and bovine feces, bovine rumen, anaerobic food waste digesters, and natural environments that included both aerobic and anaerobic municipal wastewater treatment plant samples). These clusters were primarily molded by ecology and not geography, similar to results previously reported (Fondi et al., 2016). The bovine and human feces β-lactamase clusters shared relatively high levels of connectivity, suggesting exchange of ARGs between these environments (via transfer of bacteria or horizontal gene transfer). Several previously published reports describing β-lactamase genes in bacteria from food-producing animals (Hasman et al., 2005; Cloeckaert et al., 2007; Smet et al., 2008; Costa et al., 2009), companion animals (Costa et al., 2004; Carattoli et al., 2005; Moreno et al., 2008), and wild animals (Costa et al., 2006; Poeta et al., 2008) with high similarity to β-lactamases from clinically associated microbiomes support this phenomenon. Four β-lactamase genes (blaEBR, cfxA, mecA HGI) previously identified in human-associated pathogenic and commensal Gram-positive and Gram-negative bacteria (Bellais et al., 2002; Avelar et al., 2003; Ferreira et al., 2007; Sommer et al., 2009; Hu et al., 2013; Ballhausen et al., 2014) in hospital and municipal wastewater (Schwartz et al., 2003; Bockelmann et al., 2008), were ubiquitous to these environments based on indicator analyses. In both human and bovine feces, cfxA homologs dominated the β-lactamase gene pool, indicating a selective advantage of this gene in the mammal feces. This gene is ubiquitous to Bacteroidetes isolates from anaerobic human-associated environments (Meggese and Abratt, 2015; Binta and Patel, 2016), which as describe above was highly copious in both the human and bovine feces analyzed in this study. The fact that cfxA, is frequently associated with the Tn4555 transposon (Agga et al., 2015), indicates its capacity to be horizontally transferred between phylogenetically related strains. This suggests that when evaluating microbial transfer between animals and humans, specific ARGs should be monitored concurrent to evaluation of zoonotic pathogens.
In contrast to the relatively strong interactions observed between the bovine and human fecal metagenomes, almost no network connectivity occurred between either of the fecal clusters and “environmental” cluster I. These findings are in accordance with metagenomics (Agga et al., 2015; Li et al., 2015; Munck et al., 2015; Fondi et al., 2016) and more classical approaches (Smillie et al., 2011; Forsberg et al., 2014), which reported specificity of ARGs to particular environments, with a high degree of intra-habitat mobility and low probability of ARG transfer between different environments. The significant differences in bacterial community composition between the “fecal” and “environmental” samples suggests that horizontal transfer of ARGs is often constrained by phylogenetic boundaries as previously indicated for the soil microbiome (Forsberg et al., 2014). Furthermore, other factors such as founder effects, ecological connectivity between donor and acceptor, fitness costs or second-order selection may limit the gene transfer of ARGs between environments (Martinez, 2012; Perry and Wright, 2013; Forsberg et al., 2014).
Specific evaluation of blaTEM, blaCTX-M, blaGES, and blaOXA (β-lactamases characteristic of certain Gram-negative pathogens) in the targeted environments, revealed high abundance of blaCTX-M in soils, and ubiquitous distribution of blaOXA genes across all of the analyzed environments. The similarity between the environmental β-lactamases and the clinically derived database sequences was generally low (50–60%), supporting the hypothesis that natural environments may be reservoirs of precursor ARGs that can eventually emerge in human pathogens (D’Costa et al., 2006). Conversely, these low levels of identity also indicate that pathogen-associated β-lactamases are not profuse in natural environments. Albeit the resistome hypothesis (D’Costa et al., 2007), the low connectivity observed between natural β-lactamase clusters and fecal clusters in the network analysis suggest that these events occur at relatively low frequencies. This is supported by several studies (Riesenfeld et al., 2004; Forsberg et al., 2014; Hatosy and Martiny, 2015), which also suggest limited mobility of ARGs between environmental and human-associated microbiomes. This analysis targeted a limited scope of clinically important β-lactamases, and undoubtedly, analysis of other important genes such as blaSHV need be performed in the future.
The findings presented here provide a holistic overview of β-lactamase homologs across a diverse array of environments. Nonetheless, the results should be taken with caution because screening of metagenomic data only revealed the most dominant β-lactamases in the targeted environments, and therefore evidence of gene transfer for less abundant genes will inevitably be overlooked using this approach. Furthermore, data from public metagenome projects is often highly heterogeneous, and although all of the targeted metagenomes analyzed in this study used Illumina technology, each project employed different methodologies and protocols, which may influence the obtained results. For instance, the marine metagenomes from the Sydney Harbor project (Supplementary Table S1) were sampled at 0.3 m depth, whereas the metagenomes from the Amazon continuum project were obtained from 3.63 to 4.47 m depth, at different distances from the coast. This may cause substantial differences in acquired resistomes, based on previous studies that have shown differences in ARG profiles at different water depths (Mudryk and Skórczewski, 2000; Di Cesare et al., 2015). In addition, agronomic practices can significantly affect the resistomes of agricultural soils as implied in a recent study that demonstrated that manure application stimulates shifts in soil bacterial community and ARG composition (Agga et al., 2015). In the case of feces metagenomes, the health condition and medicated status of the sampled individuals can clearly affect both the microbiological and genomic content. While ubiquity, power transformations and outlier removal partially reduces these effects, the lack of metadata associated with metagenomes can substantially complicate data interpretation.
Collectively, this study clearly indicates that both ecology and phylogeny are major driving forces in determining β-lactamase distribution in the environment, and that natural ecosystems are important reservoirs of β-lactamase gene homologs. Despite the potential exchange of β-lactamases within individual environments and between natural ecosystems and wastewater facilities, our results suggests that the exchange of β-lactamases between natural environments and human and bovine fecal microbiomes occurs at low frequencies, and therefore these events are difficult to detect. Finally, our results supports previous studies that suggest limited mobility of ARGs from natural to human-associated environments. Notwithstanding the above, metagenomics analyses only pinpoint profuse genes and we cannot rule out the notion that rare events may facilitate the transfer of novel β-lactamases from natural environments to human/animal microbiomes, thereby contributing to the propagation of AR. Future studies should focus on pinpointing these events and assessing the epidemiological risks associated with them.
Statements
Author contributions
JG drafted the main manuscript and performed the data analysis. EJ and EC were responsible for guiding and supporting the experiment and manuscript revisions.
Funding
This study was funded by the Israel Ministry of Agriculture Chief Scientist (grant # 821-0000-13). The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
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.2019.00146/full#supplementary-material
Footnotes
1.^http://metagenomics.anl.gov/
2.^https://www.ebi.ac.uk/metagenomics/
3.^http://app.agri.gov.il/eddie/EX-B_Database/EX-B_Database.fasta
References
1
AggaG.ArthurT.DursoL.HarhayD.SchmidtJ. W. (2015). Antimicrobial-resistant bacterial populations and antimicrobial resistance genes obtained from environments impacted by livestock and municipal waste.PLoS One10:e0132586. 10.1371/journal.pone.0132586
2
AllenH. K.MoeL. A.RodbumrerJ.GaarderA.HandelsmanJ. (2009). Functional metagenomics reveals diverse beta-lactamases in a remote Alaskan soil.ISME J.3243–251. 10.1038/ismej.2008.86
3
AminovR. I. (2009). The role of antibiotics and antibiotic resistance in nature.Environ. Microbiol.112970–2988. 10.1111/j.1462-2920.2009.01972.x
4
AmosG. C.ZhangL.HawkeyP. M.GazeW. H.WellingtonE. M. (2014). Functional metagenomic analysis reveals rivers are a reservoir for diverse antibiotic resistance genes.Vet. Microbiol.171441–447. 10.1016/j.vetmic.2014.02.017
5
AvelarK. E.OtsukiK.VicenteA. C.VieiraJ. M.de PaulaG. R.DominguesR. M.et al (2003). Presence of the cfxA gene in Bacteroides distasonis.Res. Microbiol.154369–374. 10.1016/S0923-2508(03)00093-7
6
BallhausenB.KriegeskorteA.SchleimerN.PetersG.BeckerK. (2014). The mecA homolog mecC confers resistance against β-lactams in Staphylococcus aureus irrespective of the genetic strain background.Antimicrob. Agents Chemother.583791–3798. 10.1128/AAC.02731-13
7
BaltzR. H. (2008). Renaissance in antibacterial discovery from actinomycetes.Curr. Opin. Pharmacol.8557–563. 10.1016/j.coph.2008.04.008
8
BellaisS.GirlichD.KarimA.NordmannP. (2002). EBR-1 a novel Ambler subclass B1 beta-lactamase from Empedobacter brevis.Antimicrob. Agents Chemother.463223–3227. 10.1128/AAC.46.10.3223-3227.2002
9
BerglundB. (2015). Environmental dissemination of antibiotic resistance genes and correlation to anthropogenic contamination with antibiotics.Infect. Ecol. Epidemiol.528564. 10.3402/iee.v5.28564
10
BhullarK.WaglechnerN.PawlowskiA.KotevaK.BanksE. D.JohnstonM. D.et al (2012). Antibiotic resistance is prevalent in an isolated cave microbiome.PLoS One7:e34953. 10.1371/journal.pone.0034953
11
BintaB.PatelM. (2016). Detection of cfxA2, cfxA3, and cfxA6 genes in beta-lactamase producing oral anaerobes.J. Appl. Oral Sci.24142–147. 10.1590/1678-775720150469
12
BockelmannU.DorriesH. H.Ayuso-GabellaM. N.de MarcayM. S.TandoiV.LevantesiC.et al (2008). Quantitative PCR monitoring of antibiotic resistance genes and bacterial pathogens in three European artificial groundwater recharge systems.Appl. Environ. Microbiol.75154–163. 10.1128/AEM.01649-08
13
CarattoliA.LovariS.FrancoA.CordaroG.MatteoP. D.BattistiA. (2005). Extended-spectrum β-lactamases in Escherichia coli isolated from dogs and cats in Rome, Italy, from 2001 to 2003.Antimicrob. Agents Chemother.49833–835. 10.1128/AAC.49.2.833-835.2005
14
CloeckaertA.PraudK.DoubletB.BertiniA.CarattoliA.ButayeP.et al (2007). Dissemination of an extended-spectrum-β-lactamase blaTEM-52 gene-carrying IncI1 plasmid in various Salmonella enterica serovars isolated from poultry and humans in Belgium and France.Antimicrob. Agents Chemother.511872–1875. 10.1128/AAC.01514-06
15
CostaD.PoetaP.BrinasL.SaenzL.RodriguesJ.TorresC. (2004). Detection of CTX-M-1 and TEM-52 β-lactamases in Escherichia coli strains from healthy pets in Portugal.J. Antimicrob. Chemother.54960–961. 10.1093/jac/dkh444
16
CostaD.PoetaP.SaenzY.VinueL.Rojo-BezaresB.JouniA.et al (2006). Detection of Escherichia coli harbouring extended-spectrum β-lactamases of the CTX-M, TEm and SHV clases in faecal samples of wild animals in Portugal.J. Antimicrob. Chemother.581311–1312. 10.1093/jac/dkl415
17
CostaD.VinuéL.PoetaP.CoelhoA. C.MatosM.SaénzY.et al (2009). Prevalence of extended-spectrum beta-lactamase producing Escherichia coli isolates in faecal samples of broilers.Vet. Microbiol.138339–344. 10.1016/j.vetmic.2009.03.029
18
Cristóbal-AzkarateJ.DunnJ. C.DayJ. M. W.Amábile-CuevasC. F. (2014). Resistance to antibiotics of clinical relevance in the fecal microbiota of mexican wildlife.PLoS One9:e107719. 10.1371/journal.pone.0107719
19
CzekalskiN.SigdelR.BirtelJ.MatthewsB.BürgmannH. (2015). Does human activity impact the natural antibiotic resistance background? Abundance of antibiotic resistance genes in 21 Swiss lakes.Environ. Int.8145–55. 10.1016/j.envint.2015.04.005
20
DaviesJ. (1994). Inactivation of antibiotics and the dissemination of resistance genes.Science264375–382. 10.1126/science.8153624
21
DaviesJ.DaviesD. (2010). Origins and evolution of antibiotic resistance.Microbiol. Mol. Biol. Rev.74417–433. 10.1128/MMBR.00016-10
22
D’CostaV.McGrannK.HughesD.WrightG. (2006). Sampling the antibiotic resistome.Science311374–377. 10.1126/science.1120800
23
D’CostaV. M.GriffithsE.WrightG. D. (2007). Expanding the soil antibiotic resistome: exploring environmental diversity.Curr. Opin. Microbiol.10481–489. 10.1016/j.mib.2007.08.009
24
D’CostaV. M.KingC. E.KalanL.MorarM.SungW. W.SchwarzC.et al (2011). Antibiotic resistance is ancient.Nature477457–461. 10.1038/nature10388
25
DelfosseV.GirardE.BirckC.DelmarcelleM.DelarueM.PochO.et al (2009). Structure of the archaeal Pab87 peptidase reveals a novel self-compartmentalizing protease family.PLoS One4:e4712. 10.1371/journal.pone.0004712
26
Di CesareA.EckertE. M.TeruggiA.FontanetoD.BertoniR.CallieriC.et al (2015). Constitutive presence of antibiotic resistance genes within the bacterial community of a large subalpine lake.Mol. Ecol.243888–3900. 10.1111/mec.13293
27
DonatoJ. J.MoeL. A.ConverseB. J.SmartK. D.BerkleinF. C.McManusP. S.et al (2010). Metagenomic analysis of apple orchard soil reveals antibiotic resistance genes encoding predicted bifunctional proteins.Appl. Environ. Microbiol.764396–4401. 10.1128/AEM.01763-09
28
El SalabiA.WalshT.ChouchaniC. (2013). Extended spectrum β-lactamases, carbapenemases, and mobile genetic elements responsible for antibiotic resistance in Gram-negative bacteria.Crit. Rev. Microbiol.39113–122. 10.3109/1040841X.2012.691870
29
European Centre for Disease Prevention and Control [ECDC] (2009). Surveillance of Antimicrobial consumption in Europe,2012. Stockholm: ECDC, 10.2900/32937
30
European Surveillance of Antimicrobial Consumption [ESAC] (2009). In ESAC yearbook 2009. Document prepared by ESAC management team, the ESAC scientific advisory board and the ESAC national networks.Stockholm: ESAC.
31
FerreiraL.AvelarK.VieiraJ.de PaulaG.ColomboA.DominguesR.et al (2007). Association between the cfxA gene and transposon Tn4555 in Bacteroides distasonis strains and other Bacteroides species.Curr. Microbiol.54348–353. 10.1007/s00284-006-0411-0
32
FondiM.KarkmanA.TamminenM. V.BosiE.VirtaM.FaniR.et al (2016). Every gene is everywhere but the environment selects: global geolocalization of gene sharing in environmental samples through network analysis.Genome Biol. Evol.81388–1400. 10.1093/gbe/evw077
33
ForsbergK. J.PatelS.GibsonM.LauberC. L.KnightR.FiererN.et al (2014). Bacterial phylogeny structures soil resistomes across habitats.Nature509612–616. 10.1038/nature13377
34
ForsbergK. J.ReyesA.WangB.SelleckE.SommerM. O.DantasG. (2012). The shared antibiotic resistome of soil bacteria and human pathogens.Science3371107–1111. 10.1126/science.1220761
35
GarauG.Di GuilmiA. M.HallB. G. (2005). Structure-based phylogeny of the metallo-β-lactamases.Antimicrob. Agents Chemother.492778–2784. 10.1128/AAC.49.7.2778-2784.2005
36
GaticaJ.TripathiV.GreenS.ManaiaC. M.BerendonkT.CacaceD.et al (2016). High throughput analysis of integron gene cassettes in wastewater environments.Environ. Sci. Technol.5011825–11836. 10.1021/acs.est.6b03188
37
GijónD.CuriaoT.BaqueroF.CoqueT. M.CantónR. (2012). Fecal carriage of carbapenemase-producing Enterobacteriaceae: a hidden reservoir in hospitalized and nonhospitalized patients.J. Clin. Microbiol.501558–1563. 10.1128/JCM.00020-12
38
GishW.StatesD. J. (1993). Identification of protein coding regions by database similarity search.Nat. Genet.3266–272. 10.1038/ng0393-266
39
HallB. G.BarlowM. (2004). Evolution of the serine beta-lactamases: past present and future.Drug Resist. Updat.7111–123. 10.1016/j.drup.2004.02.003
40
HasmanH.MeviusD.VeldmanK.OlesenI.AarestrupF. M. (2005). β-Lactamases among extended-spectrum β-lactamase (ESBL)-resistant Salmonella from poultry, poultry products and human patients in the Netherlands.J. Antimicrob. Chemother.56115–121. 10.1093/jac/dki190
41
HatosyS. M.MartinyA. C. (2015). The ocean as a global reservoir of antibiotic resistance genes.Appl. Environ. Microbiol.817593–7599. 10.1128/AEM.00736-15
42
HuY.YangX.QinJ.LuN.ChengG.WuN.et al (2013). Metagenome-wide analysis of antibiotic resistance genes in a large cohort of human gut microbiota.Nat. Commun.4:2151. 10.1038/ncomms3151
43
KellyJ.MoewsP.KnoxJ.FrêreJ. M.GhuysenJ. M. (1982). Penicillin target enzyme and the antibiotic binding site.Science218479–481. 10.1126/science.7123246
44
KiveläM.Arnaud-HaondS.SaramäkiJ. (2014). EDENetworks: a user-friendly software to build and analyze networks in biogeography ecology and population genetics.Mol. Ecol. Resour.15117–122. 10.1111/1755-0998.12290
45
LiB.YangY.MaL.JuF.GuoF.TiedjeJ. M.et al (2015). Metagenomic and network analysis reveal wide distribution and co-occurrence of environmental antibiotic resistance genes.ISME J.92490–2502. 10.1038/ismej.2015.59
46
MartinezJ. L. (2012). Bottlenecks in the transferability of antibiotic resistance from natural ecosystems to human bacterial pathogens.Front. Microbiol.2:265. 10.3389/fmicb.2011.00265
47
MassovaI.MobasheryS. (1998). Kinship and diversification of bacterial penicillin-binding proteins and β-lactamases.Antimicrobs Agents Chemother.421–17. 10.1128/AAC.42.1.1
48
McCuneB.MeffordM. (2011). PC-ORD. Multivariate Analysis of Ecological Data. Version6.Gleneden Beach, OR: User’s booklet. MjM Software.
49
MedeirosA. A. (1997). Evolution and dissemination of β-lactamases accelerated by generations of β-lactam antibiotics. In Recent increases in resistance: mechanisms and organisms.Clin. Infect. Dis.24(Suppl. 1), S19–S45. 10.1093/clinids/24.Supplement_1.S19
50
MeggeseR.AbrattV. (2015). The occurrence of antibiotic resistance genes in drug resistant Bacteroides fragilis isolates from Groote Schuur Hospital, South Africa.Anaerobe321–6. 10.1016/j.anaerobe.2014.11.003
51
MorenoA.BelloH.GuggianaD.DominguezM.GonzalezG. (2008). Extended-spectrum β-lactamases belonging to CTX-M group produced by Escherichia coli strains isolated from companion animals treated with enrofloxacin.Vet. Microbiol.129203–208. 10.1016/j.vetmic.2007.11.011
52
MudrykZ. J.SkórczewskiP. (2000). Occurrence and activity of lipolytic bacterioneuston and bacterioplankton in the estuarine Lake Gardno.Estuar. Coast. Shelf. Sci.51763–772. 10.1006/ecss.2000.0719
53
MunckC.AlbertsenM.TelkeA.EllabaanM.NielsenP. H.SommerM. O. (2015). Limited dissemination of the wastewater treatment plant core resistome.Nat. Commun.6:8452. 10.1038/ncomms9452
54
NesmeJ.CecillonS.DelmontT.MonierJ.VogelT.SimonetP. (2014). Large-scale metagenomic-based study of antibiotic resistance in the environment.Curr. Biol.241096–1100. 10.1016/j.cub.2014.03.036
55
PallecchiL.LucchettiC.BartoloniA.BartalesiF.MantellaA.GamboaH.et al (2007). Population structure and resistance genes in antibiotic-resistant bacteria from a remote community with minimal antibiotic exposure.Antimicrob. Agents Chemother.511179–1184. 10.1128/AAC.01101-06
56
PatersonD. L.BonomoR. A. (2005). Extended-spectrum beta-lactamases: a clinical update.Clin. Microbiol. Rev.18657–686. 10.1128/CMR.18.4.657-686.2005
57
PerryJ. A.WrightG. D. (2013). The antibiotic resistance mobilome: searching for the between environment and clinic.Front. Microbiol.4:138. 10.3389/fmicb.2013.00138
58
PoetaP.RadhouaniH.IgrejasG.GonçalvesA.CarvalhoC.RodriguesJ.et al (2008). Seagulls of the Berlengas Natural Reserve of Portugal as carriers of fecal Escherichia coli harboring CTX-M and TEM extended-spectrum beta-lactamases.Appl. Environ. Microb.747439–7441. 10.1128/AEM.00949-08
59
PopowskaM.RzeczyckaM.MiernikA.Krawczyk-BalskaA.WalshF.DuffyB. (2012). Influence of soil use on prevalence of tetracycline, streptomycin, and erythromycin resistance and associated resistance genes.Antimicrob. Agents Chemother.561434–1443. 10.1128/AAC.05766-11
60
RamphalR.AmbroseP. (2006). Extended-spectrum beta-lactamases and clinical outcomes: current data.Clin. Infect. Dis.42S164–S172. 10.1086/500663
61
RiceL. (2001). Evolution and clinical importance of extended-spectrum β-lactamases.Chest119391S–396S. 10.1378/chest.119.2_suppl.391S
62
RiesenfeldC. S.GoodmanR. M.HandelsmanJ. O. (2004). Uncultured soil bacteria are a reservoir of new antibiotic resistance genes.Environ. Microbiol.6981–989. 10.1111/j.1462-2920.2004.00664.x
63
SchwartzT.KohnenW.JansenB.ObstU. (2003). Detection of antibiotic-resistance bacteria and their resistance genes in wastewater, surface water and drinking water biofilms.FEMS Microbiol. Ecol.43325–335. 10.1111/j.1574-6941.2003.tb01073.x
64
ShahidM. (2014). Prevalence of CTX M extended-spectrum beta-lactamases in clinical gram-negative bacteria.Bahrain Med. Bull.36228–231. 10.12816/0008130
65
ShannonP. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks.Genome Res.132498–2504. 10.1101/gr.1239303
66
SmetA.MartelA.PersoonsD.DewulfJ.HeyndrickxM.CatryB.et al (2008). Diversity of extended-spectrum β-lactamases and class C β-lactamases among cloacal Escherichia coli in Belgian broiler farms.Antimicrob. Agents Chemother.521238–1243. 10.1128/AAC.01285-07
67
SmillieC. S.SmithM. B.FriedmanJ.CorderoO. X.DavidL. A.AlmE. J. (2011). Ecology drives a global network of gene exchange connecting the human microbiome.Nature480241–244. 10.1038/nature10571
68
SommerM. O.DantasG.ChurchG. M. (2009). Functional characterization of the antibiotic resistance reservoir in the human microflora.Science3251128–1131. 10.1126/science.1176950
69
SuJ. Q.WeiB.XuC.QiaoM.ZhuY. G. (2014). Functional metagenomic characterization of antibiotic resistance genes in agricultural soils from China.Environ. Int.659–15. 10.1016/j.envint.2013.12.010
70
SullivanR.SchausD. (2015). Extended spectrum beta-lactamases: a minireview of clinical relevant groups.J. Med. Microbiol. Diagn.4:203. 10.4172/2161-0703.1000203
71
Torres-CortésG.MillánV.Ramírez-SaadH. C.Nisa-MartínezR.ToroN.Martínez-AbarcaF. (2011). Characterization of novel antibiotic resistance genes identified by functional metagenomics on soil samples.Environ. Microbiol.131101–1114. 10.1111/j.1462-2920.2010.02422.x
72
Van BoeckelT. P.GandraS.AshokA.CaudronQ.GrenfellB. T.LevinS. A.et al (2014). Global antibiotic consumption 2000 to 2010: an analysis of national pharmaceutical sales data.Lancet Infect. Dis.14742–750. 10.1016/S1473-3099(14)70780-7
73
VersluisD.D’AndreaM. M.GarciaJ. R.LeimenaM. M.HugenholtzF.ZhangJ.et al (2015). Mining microbial metatranscriptomes for expression of antibiotic resistance genes under natural conditions.Sci. Rep.5:11981. 10.1038/srep11981
74
VilledieuA.Diaz-TorresM. L.HuntN.McNabR.SprattD. A.WilsonM.et al (2003). Prevalence of tetracycline resistance genes in oral bacteria.Antimicrob. Agents Chemother.47878–882. 10.1128/AAC.47.3.878-882.2003
75
WagnerU.PetersenE.SchwabH.KratkyC. (2002). EstB from Burkholderia gladioli: a novel esterase with a β-lactamase fold reveals steric factors to discriminate between esterolytic and β-lactam cleaving activity.Protein Sci.11467–478. 10.1110/ps.33002
76
WalshF.DuffyB. (2013). The culturable soil antibiotic resistome: a community of multi-drug resistant bacteria.PLoS One12:e65567. 10.1371/journal.pone.0065567
77
YangJ.WangC.ShuC.LiuL.GengJ.HuS.et al (2013). Marine sediment bacteria harbor antibiotic resistance genes highly similar to those found in human pathogens.Microb. Ecol.65975–981. 10.1007/s00248-013-0187-2
78
YangY.JiangX.ChaiB.MaL.LiB.ZhangA.et al (2016). ARGs-OAP: online analysis pipeline for antibiotic resistance genes detection from metagenomics data using an integrated structured ARG-database.Bioinformatics322346–2351. 10.1093/bioinformatics/btw136
Summary
Keywords
antibiotic resistance, β-lactamases, network analysis, metagenomics, environment
Citation
Gatica J, Jurkevitch E and Cytryn E (2019) Comparative Metagenomics and Network Analyses Provide Novel Insights Into the Scope and Distribution of β-Lactamase Homologs in the Environment. Front. Microbiol. 10:146. doi: 10.3389/fmicb.2019.00146
Received
06 July 2018
Accepted
21 January 2019
Published
11 February 2019
Volume
10 - 2019
Edited by
Patrick Rik Butaye, Ross University School of Veterinary Medicine, Saint Kitts and Nevis
Reviewed by
Takehiko Kenzaka, Osaka Ohtani University, Japan; Katharina Schaufler, University of Greifswald, Germany
Updates
Copyright
© 2019 Gatica, Jurkevitch and Cytryn.
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: Eddie Cytryn, eddie@volcani.agri.gov.il
This article was submitted to Antimicrobials, Resistance and Chemotherapy, a section of the journal Frontiers in Microbiology
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.