A sampling and metagenomic sequencing-based methodology for monitoring antimicrobial resistance in swine herds

Objectives

Reliable methods for monitoring antimicrobial resistance (AMR) in livestock and other reservoirs are essential to understand the trends, transmission and importance of agricultural resistance. Quantification of AMR is mostly done using culture-based techniques, but metagenomic read mapping shows promise for quantitative resistance monitoring.

Methods

We evaluated the ability of: (i) MIC determination for Escherichia coli; (ii) cfu counting of E. coli; (iii) cfu counting of aerobic bacteria; and (iv) metagenomic shotgun sequencing to predict expected tetracycline resistance based on known antimicrobial consumption in 10 Danish integrated slaughter pig herds. In addition, we evaluated whether fresh or manure floor samples constitute suitable proxies for intestinal sampling, using cfu counting, qPCR and metagenomic shotgun sequencing.

Results

Metagenomic read-mapping outperformed cultivation-based techniques in terms of predicting expected tetracycline resistance based on antimicrobial consumption. Our metagenomic approach had sufficient resolution to detect antimicrobial-induced changes to individual resistance gene abundances. Pen floor manure samples were found to represent rectal samples well when analysed using metagenomics, as they contain the same DNA with the exception of a few contaminating taxa that proliferate in the extraintestinal environment.

Conclusions

We present a workflow, from sampling to interpretation, showing how resistance monitoring can be carried out in swine herds using a metagenomic approach. We propose metagenomic sequencing should be part of routine livestock resistance monitoring programmes and potentially of integrated One Health monitoring in all reservoirs.

Introduction

Antimicrobial resistance (AMR) is considered one of the greatest threats to human health.¹ The resistance epidemic has been attributed to the use of antimicrobials in clinical settings, and in livestock.² In addition to the abundance of resistance in animal reservoirs, animals are also important in the early evolution and emergence of novel resistance genes.³

Many research studies have shown that the use of antimicrobials in livestock will lead to increased occurrence of AMR.^4–6 In addition, research shows that reducing the use of antimicrobials can decrease the occurrence of resistance.^7,8 In order to identify areas of priority and document effects of interventions, a robust monitoring system is essential. The first integrated routine monitoring of AMR was established in Denmark in 1995.⁹ Since its inception, other nations have developed similar monitoring programmes that embrace the One Health approach.^10,11

Current monitoring efforts are mainly based on culturing indicator bacteria and determining their resistance phenotypically.¹² The procedures for sampling and analysis are time consuming, expensive and only target a limited number of the bacteria and resistance mechanisms found in the different reservoirs.

Traditionally, the proportion of resistant bacteria has been assessed through cfu counting, with and without a selective agent, and used as a measure of resistance.⁶ Alternatively, phenotypic MIC determination provides high accuracy, reproducibility and the isolate-specific context to monitor multiresistance. Both methods are, however, blind to the cause of resistance and the epidemiology of the responsible genes. To alleviate this, isolated bacteria can be screened with PCR using primers targeting the suspected causal resistance genes.^13,14 More recently, WGS of indicator organisms has become economically feasible and pipelines for in silico resistance determination have become accurate enough to compete with phenotypic resistance determination and PCR.¹⁵ For monitoring purposes, however, a prohibitively large number of isolates needs to be analysed to obtain reliable prevalence estimates for resistance.

Indicator-organism-centric approaches only provide insight into a minor fraction of naturally occurring bacteria in which the bulk of resistance genes may be present.¹⁶ Metagenomic read mapping has recently been used to quantify resistance genes directly in multiple reservoirs, including waste water, humans, livestock, drinking water, activated sludge and aircraft septic tanks.^17–21 The approach has furthermore proved useful in quantifying the selection pressure induced by ciprofloxacin in infants.²² Using metagenomics, one can quantify thousands of targets in a sample without requiring a priori knowledge of which genes or bacteria are present. A single dataset of DNA sequences can therefore be used to quantify all known resistance genes.

The aim of our study is to compare methods for measuring the occurrence of AMR in swine herds. In particular, we focus on sampling and analytic procedures, and we test whether pen floor sampling provides a valid approximation to sampling individual pigs. We compare metagenomics with traditional culture-based methods for quantifying AMR in swine herds. The study was carried out in 10 herds.

Materials and methods

Choice of sample herds

In Denmark, detailed records are kept on all antimicrobial prescriptions for livestock herds. The farmer obtains drug prescriptions from veterinarians and the information is recorded in the VETSTAT database.²³ Another database, the Danish Central Husbandry Register (CHR), collects mandatory information from all livestock herds including ownership, animal counts, animal species and age groups (e.g. sows/piglets, weaners, finishers). Using these databases, we generated a list of herds fitting our study population criteria: conventional Danish integrated pig herds with more than 500 sows and producing at least 5000 slaughter pigs annually. For these herds, we calculated the defined animal daily doses (ADD) for antimicrobial drug groups and adjusted for the CHR-based herd sizes.²⁴ To include herds covering a wide spectrum of antimicrobial consumption, we invited herd owners from the top and bottom 10% tetracycline consumption quantiles to participate in the study, until five herds within each quantile had accepted. For information on which active compounds make up each drug class category, see Table S1 (available as Supplementary data at JAC Online).

Sampling and pooling

For an overview of sampling, pooling and downstream analysis, see Figure S1. Briefly, from each of the 10 herds, herd-level floor (HL-F) samples were obtained from 30 pens and paired pen-level floor (PL-F) and pig (PL-P) samples were obtained from four selected pens. Within each pen, the sampled pigs were randomly chosen. For details, see the Supplementary data.

Cultivation of faecal bacteria

For each pooled sample, 1 g of faeces was suspended in 9 mL of isotonic saline and serially diluted. Aliquots of 100 µL dilutions were plated onto selective and non-selective LB and MacConkey agar plates to quantify aerobic bacteria and Escherichia coli, respectively. Selective plates contained 8 mg/L tetracycline (T3383 tetracycline hydrochloride, Sigma–Aldrich) in LB and 16 mg/L tetracycline or ampicillin (A9393 Ampicillin, Sigma–Aldrich) in MacConkey. All assays were performed in triplicate and plates were incubated at 37°C overnight. For each triplicate set, a weighted average resistance proportion was calculated. For details, see the Supplementary data.

MIC

We determined MIC values for a subset of E. coli from the cultivation experiments using the previously described DKMVN4 panel.²⁵ For details, see the Supplementary data.

DNA extraction

For DNA extraction, a modified QIAamp Fast DNA Stool Mini Kit protocol was employed (Qiagen, Valencia, CA, USA). Cryotubes with pooled faeces were gently thawed on ice. Prior to the protocol, 0.2 g of sample was mixed with 1 mL of InhibitEX buffer in a Lysing Matrix A tube (MP Biomedicals). Samples were treated with a TissueLyser (3 × 30 s, 30 Hz) and were chilled on ice between repetitions. Following bead beating, samples were heated to 95°C for 7 min and centrifuged to eliminate stool particles. During the extraction protocol, volumes of sample, proteinase K, buffer AL and ethanol were doubled. DNA was eluted in 100 µL of elution buffer.

Shotgun DNA sequencing

PCR-free DNA libraries were generated and sequenced on the HiSeq2500 (Illumina) to generate roughly 7 gigabases of paired-end reads per sample, enough to get 20× coverage of bacteria with 1% abundance.²⁶ Cutadapt was used to trim the reads to a mean Phred score of 20.²⁷ For details, see the Supplementary data.

Read mapping (resistance)

Resistance was quantified using the ResFinder database, a collection of 2130 annotated resistance genes.¹⁵ Paired-end mapping was done with a BWA-mem-based script as described elsewhere.^20,28 To avoid favouring poorer-quality samples (those that on average get shorter reads post-trimming), we required each read to align with a static 50 bp instead of the default setting that is relative to the trimmed read length (80%).

Because the ResFinder database contains many highly identical sequences, unspecific mapping occurs when reads map to identical parts of homologous gene variants. Read counts from variants of the same gene were aggregated to gene levels according to common gene names.²⁰ For the variant-to-common name conversion table and the final abundance matrix, see Tables S2 and S3, respectively.

Read mapping (bacteria)

For the bacterial count data, we employed the EBI Metagenomics resource.^29,30 Briefly, the pipeline identifies DNA reads matching the bacterial 16S gene. It then uses QIIME to map those to the Greengenes 16S database and to identify the lowest common ancestor of each 16S sequence.^31,32 The taxonomic abundance matrices can be downloaded from EBI Metagenomics and are associated with the accession number (ERP013942).

Metagenomic data analysis

A resistance gene count matrix was constructed from the DNA read mapping output, with genes (2130) and samples (20) constituting rows and columns, respectively (Table S3). Values >0 were divided by 2 so each count reflected a single mapping fragment. Read counts were aggregated according to the common gene names specified in Table S2, resulting in the final abundance matrix. To avoid inappropriate statistical models and losing power by rarefying, metagenomic data such as RNAseq data, should be analysed using zero-inflated negative binomial models as suggested by McMurdie and Holmes.³³ Therefore, gene-level differential abundance analysis was carried out using the R-package DESeq2.³⁴ For details on differential abundance analysis, see the Supplementary data.

Expected resistance

The expected tetracycline resistance was derived by ranking the farms from 1 to 10 in terms of VETSTAT-derived, herd-size-adjusted tetracycline consumption 1 year prior to sampling. Ranks based on resistance measurements can thus be compared with the expected resistance, based on the assumption that tetracycline resistance occurrence follows tetracycline consumption.

qPCR

Primers and probes have been presented previously by Schmidt et al.³⁵ The inclusion criteria for the AMR genes were: (i) the use of the antimicrobial class linked to the AMR gene in question in the Danish pig production; (ii) the occurrence of the gene in a wide bacterial population; and (iii) the possibility of designing a qPCR assay for the chosen genes utilizing the same temperature profile.³⁵ Quantification of the four AMR genes [erm(B), sul1, tet(B) and tet(W)] was done using the high-capacity qPCR chip ‘192.24 Gene Expression IFC’ (Fluidigm^®) using two technical replicates as previously described by Clasen et al.³⁶ The efficiency of the primers was determined using standard curves and obtained results were normalized with 16S ribosomal DNA (reference gene). Negative ΔΔCq values, indicating no observed gene, were treated as missing values. Subsequently, ΔΔCq values were unlogged (2^−ΔΔCq) to directly compare with alternative resistance measurements.

Comparison of metagenomics, phenotypic methods and expected resistance

The four resistance estimates, as well as the expected tetracycline resistance were correlated to each other (Spearman rank correlation). Using a rank-based statistic, we avoid making assumptions about the shape of any relationship. Since we found that floor and rectal samples were comparable (see the Results and discussion section), both pig (PL-P) and floor (PL-F/HL-F) isolates were used for the MIC estimate, as too few isolates were obtained from the individual HL-F samples. The remaining methods were based on the HL-F samples alone. All herd-level tetracycline resistance measurements can be found in Table S4. The bootstrapped hierarchical dendrogram was produced with the R package ‘pvclust’, using 10 000 bootstrap iterations (with default other settings), and shows node certainty percent based on the approximately unbiased P value (AUP).³⁷

Accession numbers

The DNA sequences (reads) from the 20 metagenomic samples are deposited in the European Nucleotide Archive (ENA) under the project accession number ERP013942.

Visualization

The R package RcolorBrewer was used to generate the colour palettes used in several figures. The package is based on work by Cynthia A. Brewer (http://www.ColorBrewer.org).

Results and discussion

Ten Danish swine herds with historical high and low antimicrobial consumption were enrolled in the study. To make sure their usage was also varied in the period leading up to the sampling, we recalculated their consumption for a 1 year window prior to the individual sampling dates. As expected from our selection, a large variation in overall drug use was observed between the herds (Figure 1a). This was also the case for most individual drug classes that were not considered in the initial selection of the study herds (Figure 1b and Table S5). By dosage, tetracycline was the most used drug class, followed by β-lactams and macrolides. This pattern among the study herds is consistent with Denmark's antimicrobial usage in all swine the same year.³⁸

Figure 1.

Antimicrobial consumption in the study herds. Information about antimicrobial drug usage, prescribed in a 1 year window prior to sampling, was extracted from the VETSTAT database. For each combination of herd and drug type, ADD were calculated. (a) Stacked bar chart showing total antimicrobial consumption per herd (ADD). (b) Multiple bar charts showing ADD per herd stratified by antimicrobial group.

Comparison of sampling strategies

The pig intestine is an important source of zoonotic bacteria contaminating meat during slaughter. In this study, we obtained three kinds of samples. Those were: (i) rectal pig faecal samples; (ii) undisturbed pen floor faecal samples; and (iii) pen floor manure samples. Rectal samples and manure samples are the most and least time consuming to collect, respectively. For an overview of sampling and pooling strategy, see Figure S1.

To determine whether sampling undisturbed faeces on the pen floor is an acceptable proxy for sampling recta, we compared paired PL-F and PL-P faeces pools from four pens per herd, for a total of 40 pairs. When using cfu counting, we found that tetracycline resistance correlated significantly (Spearman rank correlation, P < 0.05) between aerobic bacteria in floor and rectal samples (Figure S2a). We also found significant relationships for tetracycline-resistant (Figure S2b) and ampicillin-resistant E. coli (Figure S2c).

To confirm that pig and fresh floor faecal samples were comparable, we attempted to assess the association using a different method (qPCR) and genes representing more resistance types. We quantified erm(B), sul1, tet(B) and tet(W) in the same 40 paired sets of PL-F and PL-P samples. The erm(B) and tet(W) abundances correlated well between floor and pig samples (P < 2.2 × 10^–16) (Figure S2d and e). The genes tet(B) and sul1 were below the detection limit in several samples, yielding too few pairs for comparison. The qPCR abundances for the four genes are included in Figure S3.

While the earlier described PL-F samples were derived from fresh, undisturbed droppings, the herd-level floor (HL-F) samples originated from floor manure of unknown age, which is readily available and requires no animal handling. We sequenced DNA from the 20 herd-level samples, for a total of 690 million pairs (1380 million reads), corresponding to 139 billion bp (Table S6). This yielded an average of 7 billion bp/sample, the recommended amount to get 20× coverage on bacteria with a 1% relative abundance.²⁶ We then used those DNA sequences to quantify known resistance genes.

Principal component analysis (PCA) showed herd-wise clustering, suggesting that herd is a stronger resistance predictor than sampling location (floor versus pig) (Figure 2a). This was confirmed by hierarchical clustering based on the most common genes, where we found 7 of the 10 sample pairs clustered according to herd of origin (Figure 2b). Despite the high herd-wise similarity, we noticed a consistent pen floor shift along the second principal component. To investigate this systematic change in the resistome, we tested for differential abundance between the paired HL-P and HL-F samples. Only two resistance genes differed significantly between the sample types: tet(M) and tet(36), both overrepresented in floor manure samples (false discovery rate-adjusted P value, FDR < 0.1). tet(M) has the broadest host range of all tetracycline resistance genes, being found in both Gram-positive and Gram-negative bacteria.³⁹tet(36), on the other hand, is quite limited and primarily found in Bacteroides coprosuis, where it was first isolated from swine manure.⁴⁰ Consistent with this, the gene was not detected in any pig samples, whereas most manure samples contained the gene at low abundance. As shown by Forsberg et al.,⁴¹ different samples' resistomes are structured according to their microbial taxonomic composition. In order to check if the change in resistance was explained by a taxonomic shift, we performed differential abundance analysis on bacterial genera quantified in the same datasets. Five genera were differentially abundant (FDR < 0.1): Succinivibrio, Fusobacterium, Rummeliibacillus, Acinetobacter and Ignatzschineria (Figure S4). Like for resistance genes, all significant genera were overrepresented in the floor samples. Both Fusobacterium and Acinetobacter isolates are known to harbour tet(M) genes.^39,42Bacteroides was not differentially abundant between the sample types, but this genus also spans many species. Concordant with finding tet(36) only in floor samples, B. coprosuis-specific reads were also only identified in the floor samples (data not shown, based on EBI Metagenomics).

Figure 2.

Metagenomic associations between pig and manure floor samples. The resistance gene abundances were rlog transformed using DESeq2 prior to visualization. (a) Principal component analysis plot visualizing resistome differences, annotated according to herd of origin (colour) and sampling strategy (shape). (b) Heatmap of the resistance genes with highest average abundance across the 20 datasets. Rlog abundances are scaled to each gene (Z-score) allowing easy inter-sample comparison. Bold italic sample names indicate highest overall antimicrobial consumption.

Despite these systematic differences between the sample types, the samples clustered according to herd of origin in 9 out of 10 cases based on bacterial composition (Figure S4).

Together these findings suggest the increased resistance in floor samples is explained by an expansion of certain bacterial taxa, either through contamination or proliferation. In summary, both fresh floor droppings and manure samples represent the intestinal resistome well and can provide acceptable proxies for sampling individual pigs. Pig pen floors have earlier been sampled for risk assessment, but until recently very little was known about how well these samples represent the intestines. Multiple sampling strategies have recently been evaluated for pig herd resistance monitoring with a qPCR quantification protocol. Consistent with our results, the authors found pen floor sampling appropriate for monitoring intestinal resistance.³⁵

Comparison of resistance monitoring methods

In order to benchmark the techniques, we assumed antimicrobial use increases AMR—a frequently made observation.^5,6,18 The tetracycline ADD were used to rank the expected tetracycline resistance levels in the 10 herds relative to one another. We then used: (i) cfu counting of aerobic bacteria; (ii) cfu counting of E. coli; (iii) MIC determination for E. coli; and (iv) metagenomic read mapping to quantify the overall tetracycline resistance in the same 10 HL-F samples. We compared the methods to each other and to the expected resistance using correlation analysis and hierarchical clustering. The two phenotypic resistance measurements for E. coli (cfu- and MIC-based) correlated significantly (Spearman rank, α < 0.05) (Figure S5). This is meaningful because they both describe the same species phenotypically. See Figure S6 for all MIC panel results. The only tetracycline resistance measurement to significantly correlate with the expected resistance was the metagenomic measurement. For the hierarchical clustering of methods, the E. coli resistance estimates group together, whereas the two species-independent estimates (aerobic bacteria and metagenomics) were both closer to the expected resistance, with metagenomics being the closest (Figure S5b). Bootstrapping indicated that clustering was robust. All herd-level tetracycline resistance numbers are presented in Table S4.

The lack of negative correlations suggests that tetracycline usage increases tetracycline resistance and that this observation is independent of measurement method. However, only the metagenomic measurement significantly correlated with expected resistance. The specific correlation coefficients suggest the bulk of tetracycline resistance genes, found using metagenomic read mapping, follow tetracycline consumption closer than E. coli alone would suggest. Other studies have found significant associations between antimicrobial consumption and phenotypic resistance prevalence estimates.^5,6 We stress that our study, which has a low number of truly independent samples, does not disqualify phenotypic association studies. We simply find that the drug-resistance association is strongest with our metagenomic approach.

Gene-level associations with antimicrobial consumption

As our metagenomic measurement for tetracycline resistance was the only measurement to significantly correlate with the expected tetracycline resistance, we chose to perform regression analyses on metagenomic data to assess whether specific genes were associated with consumption of the most commonly used drug classes: tetracyclines, macrolides, narrow-spectrum penicillins and aminoglycosides.

Tetracycline use was associated with the tetracycline resistance gene tet(44) in pig samples (FDR < 0.1). Macrolide usage was associated with macrolide resistance genes erm(B) and erm(G) in floor and pig samples, respectively.

apmA, a gene providing resistance to the veterinary aminoglycoside apramycin, followed the aminoglycoside consumption, both in manure and pig samples. In the pig samples, aminoglycoside use was associated with tet(44) and two additional aminoglycoside resistance genes: ant(6)-I and str(B). In summary, associations were between resistance genes and the antibiotic drug classes they protect against. The exception was tet(44), which followed aminoglycoside consumption in pigs. tet(44) is co-located with the ant(6)-I variant ant(6)-Ib in a mobile Campylobacter pathogenicity island.⁴³ Since these both followed aminoglycoside usage, it seems likely the association is explained by aminoglycoside co-selection of tet(44). Regression analysis results are included in Table S7.

We noticed an overrepresentation of positive associations among genes with low unadjusted P values, suggesting even more genes were positively associated with drug use, but outside our resolution.

To test if this was the case, while allowing for both co- and cross-resistance, we correlated the most common resistance gene abundances in the 20 metagenomes with the drug consumption values (Figure 3a). There were many more positive correlations than negative correlations. In addition, among the positive correlations, P values <0.05 were overrepresented (8.45%) (Figure 3b). This was in contrast to the negative correlations, where low P values were underrepresented (3.21%) (Figure 3c). While no single correlation should be trusted without correcting for multiple testing, this confirms that drug use is associated with AMR, and indicates our metagenomic approach is sensitive enough to detect drug-induced changes in the swine herd resistome.

Figure 3.

Association between resistance gene abundance and drug class consumptions. Resistance gene abundances (FPKM) were Spearman correlated against consumption estimates for the most used drug classes. (a) Correlation plot where colour denotes the Spearman correlation coefficient (R) of a given correlation and the circle size is proportional to the coefficient of determination (R²). White crosses indicate P values above 0.05. (b) Histogram of P values for the positive correlation coefficients. (c) Histogram of P values for the negative correlation coefficients. (b and c) A broken line shows the P value threshold of 0.05 and, above the plots, the proportion of correlations below that threshold is indicated.

A valid concern with using metagenomic-based read mapping for measuring resistance is the inability to distinguish between zoonotic and immobile genes in an environment.⁴⁴ This drawback of metagenomics could, however, be reduced. Metagenomic assembly techniques can give contextual information on genes' locations and is an area of active development.⁴⁵ Novel approaches such as using gene co-abundance profiles across many samples, metagenomic binning or molecularly linking adjacent DNA prior to lysis can also help elucidate the taxonomic affiliation of accessory genes.^46,47 Metagenomic assembly would complement our read-based quantification well, as our proposed method is limited to quantifying resistance gene groups and is unable to identify specific alleles of interest. Rare resistance genes would not, however, get sufficient coverage for assembly at currently feasible sequencing depths, leaving an important role for read-mapping. Also, it is important to note that our chosen database ResFinder focuses on resistance genes found in culturable pathogens, meaning one might miss a large fraction of functional resistance genes. In order to determine which genes are ultimately important to human and animal health, further research is needed, especially into the host range and mobility of AMR genes, their disseminating vectors and taxonomic background. Using functional metagenomics one can detect novel resistance genes, while the aforementioned methods can help elucidate whether a discovered resistance gene is likely to transfer to potential human pathogens.^44,48 These methods should together inform which genes we need actively to screen for and add to databases. Our observational approach and more controlled experiments can then be used to determine drug resistance associations, which should inform future drug use.²²

Conclusions

We have developed an approach for measuring AMR in swine herds that can be applicable elsewhere. Key features are: (i) sampling multiple pen floors for manure, followed by pooling to a single composite sample; (ii) use of a DNA extraction kit that retains a high microbial diversity; (iii) use of PCR-free library preparation to avoid unnecessary amplification bias; (iv) paired-end sequencing for improved accuracy; and (v) quantification of DNA sequence mapping to a database of known resistance genes. Previous experimental studies have shown, using metagenomics, that therapeutic or sub-therapeutic additions of antibiotics cause increases in resistance.^49,50 To our knowledge, this is the first observational study showing that metagenomics can quantify and elucidate the effects of antimicrobial usage in livestock animals. The fact that the technology can detect these changes in natural populations makes us confident that metagenomics should be used in future AMR monitoring programmes.

Funding

This work was supported by the Danish Veterinary and Food Administration and the EU FP7 programme (grant number 613754, Effort against antimicrobial resistance, http://www.effort-against-amr.eu/). S. J. P. was supported by Carlsbergfondet (grant number 2013_01_0377). The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Transparency declarations

None to declare.

Author contributions

F. M. A., Y. A., H. V. and A. F. conceived the study. V. D. A., L. D. K., M. S. J. and H. V. designed the sampling strategy, selected study participants and performed sampling. H. M., J. C., B. E. K. and M. S. J. performed the laboratory work. P. M., H. V., F. M. A., S. J. P., Y. A. and O. L. analysed the biological and epidemiological data. P. M. prepared the manuscript. All authors reviewed and edited the manuscript.

Supplementary data

Supplementary data, including Figures S1–S6 and Tables S1–S7, are available at JAC Online (http://jac.oxfordjournals.org/).

Acknowledgements

We would like to thank the 10 farmers who participated in the study without any form of compensation.