Abundance and diversity of the fecal resistome in slaughter pigs and broilers in nine European countries

EFFORT group Haitske Graveland, Alieda van Essen, Bruno Gonzalez-Zorn, Gabriel Moyano, Pascal Sanders, Claire Chauvin, Julie David, Antonio Battisti, Andrea Caprioli, Jeroen Dewulf, Thomas Blaha, Katharina Wadepohl, Maximiliane Brandt, Dariusz Wasyl, Magdalena Skarzyñska, Magdalena Zajac, Hristo Daskalov, Helmut W Saatkamp, Katharina D.C. Stärk. Abstract Antimicrobial resistance (AMR) in bacteria and associated human morbidity and mortality is increasing. Use of antimicrobials in livestock selects for AMR that can subsequently be transferred to humans. This flow of AMR between reservoirs demands surveillance in livestock as well as in humans. As part of the EFFORT project (www.effort-against-amr.eu), we have quantified and characterized the acquired resistance gene pools (resistomes) of 181 pig and 178 poultry farms from nine European countries, generating more than 5,000 gigabases of DNA sequence, using shotgun metagenomics. We quantified acquired AMR using the ResFinder database and a database constructed for this study, consisting of AMR genes identified through screening environmental DNA. The pig and poultry resistomes were very different in abundance and composition. There was a significant country effect on the resistomes, more so in pigs than poultry. We found higher AMR loads in pigs, while poultry resistomes were more diverse. We detected several recently described, critical AMR genes, including mcr-1 and optrA, the abundance of which differed both between host species and countries. We found that the total acquired AMR level, was associated with the overall country-specific antimicrobial usage in livestock and that countries with comparable usage patterns had similar resistomes. Novel, functionally-determined AMR genes were, however, not associated with total drug use.


Introduction
Sampling to estimate the effect of random sampling 89 To study the potential effect of sampling randomness and the reproducibility of our sampling protocol, two 90 of the pig herds were chosen for triplicate sampling. These two herds were sampled three times on the 91 same day (25 samples x 3 sampling rounds), resulting in six pooled samples (2 herds x 3 sampling rounds), 92 from which the within-farm variation was assessed. A table with all samples and their metadata is included 93 as Supplementary Table 1. 94 read pairs using BWA-MEM, ensuring that the thresholds for a read pair to map as a 'proper pair' were the 121 same for all read pairs, and avoiding bad insert size estimates when few read pairs aligned. 122 The read pairs were aligned to the bacterial NCBI genomes again, and to the AMR genes present in the 123 ResFinder database (accessed Nov. 17, 2016). 20 ResFinder is a manually curated database of horizontally 124 acquired AMR genes and thus does not include intrinsic AMR genes and mutated housekeeping genes 125 providing AMR by changing the drug target. Technical duplicate read pairs were then removed using 126 'MarkDuplicates' from the Picard command line tools (v2.8.3; http://broadinstitute.github.io/picard/). 127 For each ResFinder reference sequence, we counted the number of read pairs that properly aligned with at 128 least 50 bp aligning from both the forward and reverse reads. Each read pair matching a ResFinder 129 reference was assigned to the first highest-scoring reference, as done in MGmapper. The same was done 130 7 constructed a functional resistance database (FRD) from 3,416 AMR gene variants identified in four major 151 studies, using 23 different antimicrobials for selection. 23-26 152 Briefly, in each of these studies, DNA was extracted from environmental and human fecal samples, 153 fragmented and cloned into a plasmid vector and screened for AMR functionality in E. coli cultured with 154 one of multiple antimicrobials. AMR-granting plasmid inserts were then sequenced and the responsible 155 open reading frame was identified. The protocol for the database construction can be found at 156 https://cge.cbs.dtu.dk/services/ResFinderFG/. Genes were quantified using MGmapper as was done for 157 ResFinder. Genes with more than 90% identity to ResFinder genes were removed post mapping to obtain 158 the new AMR genes without overlap with ResFinder. The resulting data was aggregated to 90% gene 159 clusters, using CD-HIT-EST, as was done for ResFinder. 22 The most frequent gene clusters remaining were 160 derived from genes selected using: trimethoprim, chloramphenicol, co-trimoxazole, cycloserine, amoxicillin, 161 gentamicin, penicillin and tetracycline.
was carried out for both pigs and poultry; combined and separately, using the vegan function 'betadisper'. 166 The same analysis was used to test whether host animal and country were significant predictors of within-167 group beta diversity dispersion. The effects of country on sample dissimilarities was determined using 168 'permutational multivariate analysis of variance using distance matrices' ('adonis2' function in vegan 169 package), separately for pig and poultry. 170

Antimicrobial use in livestock 171
Data for national livestock antimicrobial usage (AMU) was obtained from the European Medicines Agency's 172 2014 European Surveillance of Veterinary Antimicrobial Consumption (ESVAC) report and was stratified by 173 major drug family. 28 The mass of active compound sold for use in animals in 2014 was divided by the 174 Population Correction Unit (PCU) in 10 6 kg -approximating the biomass. The PCU is a unit that allows inter-175 species integration by adjusting for import/export and differences in average weight between species when 176 they are most likely to receive antimicrobial treatment. The estimate was multiplied by 1000 to obtain drug 177 mg/kg PCU livestock. The country-specific veterinary drug use can be found in Supplementary Table 5. 8 between the samples in the PCU-corrected AMU (Supplementary Table 5 The circular BC sample dendrogram was exported in Newick format using the ape package and further 201 annotated with the Interactive Tree of Life tool. 31,32 Bar-, box-and scatter plots were produced using the 202 ggplot2 R library. 33 The R library RcolorBrewer was used to generate the color palettes used. The library is 203 based on work by Cynthia A. Brewer (www.ColorBrewer.org). 204

Statistical analyses 205
All statistics was done in Microsoft R Open (MRO) 3.3.2, using the libraries and procedures detailed below. 206 Exact package versions can be found here: https://mran.revolutionanalytics.com/snapshot/2016-11-207 01/bin/windows/contrib/3.3/. For statistical tests, only the first sampling from triple-sampled herds was 208 included (see Supplementary Table 1). Unless otherwise mentioned, all statistical analyses were performed 209 on pigs and poultry separately. 210

Effect of AMU on total AMR 211
For testing the effect of total AMU on total metagenomic AMR abundance (sum of all genes), we used the 212 lme4 1.1-12 package to make linear mixed effects regression models with total livestock drug usage as the 213 independent variable, total AMR abundance (FPKM) as the dependent variable and country as a mixed 214 effect intercept, adjusting for the fact that AMR abundance observed in farms from the same country is 215 an average of 50 million (SD: 18*10^6) PE reads per pooled sample. This was similar for pig and poultry 242 samples, though the former varied more than the latter. 243 244 Acquired resistome characterization 245 The total AMR gene level varied significantly across samples, both depending on host animal and country of 246 origin. In general, pigs had a higher AMR level than poultry (Figure 1a

257
Stacked barchart of AMR abundance per type (colors) per sample (x-axis), proportional to total AMR within each sample.

258
We summed the relative abundance of AMR to the corresponding drug class level for each sample to look 259 for major trends across species and countries (Figure 1b). When considering the proportion of the total 260 resistome that each AMR type takes up, the pig samples were relatively homogenous: tetracycline AMR 261 was by far the most common, followed by macrolide AMR. Beta-lactam and aminoglycoside AMR genes 262 followed with other kinds of AMR being rare. Italian pigs had a notably larger proportion of phenicol AMR 263 compared with other countries and it seemed to be consistent across all Italian farms. A subset of Bulgarian 264 pig farms had a similar proportion of phenicol AMR. 265 Among the poultry farms, there was less consistency: both within and between countries, the relative 266 proportions of AMR per drug class varied more. Tetracycline, macrolide, beta-lactam and aminoglycoside 267 AMR made up the majority of AMR, as in pig samples, but the two latter classes had very minimal 268 contributions in a subset of herds. Sulfonamide and trimethoprim AMR were more abundant in poultry 269 samples compared with pig samples, across all countries. In many Polish poultry herds, quinolone AMR  When analyzing the two reservoirs separately, we observed clustering according to country of origin in pigs 302 (Figure 2c), while clustering was more diffuse for poultry (Figure 2d). We tested for the country effect and 303 found it to be significant in both pigs (adonis2 p<0.001) and poultry (adonis2 p<0.001). In poultry however, 304 the country effect only explained roughly a quarter of the variation, while country explained roughly half of 305 the variation in pigs (data not shown). In the pig resistome ordination, the Danish and Dutch samples 306 including cat(pC194), catP, and cat_2, were much more abundant in ltaly, compared to the other countries, 316 consistent with our inspection of AMR at class level (Figure 1). Several AMR genes known to be co-located 317 indeed co-occurred across samples. The genes in the vancomycin AMR VanA cassette were co-located in a 318 number of poultry samples. This was also true for the VanB cassette members, clustering together, but 319 separately from VanA, showing ability to distinguish variants of homologous genes. As earlier indicated, the 320 poultry samples showed less country-based clustering than pigs. An exception were the Danish poultry 321 samples. These had noticeably lower abundance of many AMR genes that were widespread in other 322 countries. 323

Core resistome 324
To determine whether specific genes were unique to each of the host animals, we examined the set of AMR 325 genes that was consistently observed within each animal species (evidence for it in 95% of samples). We 326 identified 33 core AMR genes in pigs and 49 core AMR genes in poultry, with 24 being shared between the 327 two hosts (Supplementary Figure 6). Hence, only nine AMR genes were pig-core genes without also being

Alpha diversity and richness 362
We calculated several alpha diversity indexes for each farm resistome (Figure 4) The range of AMR diversity 363 was generally much larger for poultry samples, having both lower and higher diversity than pig samples, 364 which had a tighter spread of diversity. The poultry samples had a higher estimated richness than pigs (i.e. 365 a higher number of unique AMR genes per sample). Alpha diversity indexes can be found in Supplementary 366 Table 9. 367 Interestingly, countries with high AMR richness in pigs also tended to have high AMR richness in poultry.

392
The correspondence between the two datasets was slightly stronger in pigs (Procrustes symmetric 393 correlation: 0.71) than in poultry (Procrustes symmetric correlation: 0.66). Interestingly, in pig samples we 394 saw a country effect on the strength of association between the bacteriome and the resistome. In Dutch 395 and Belgian pig herds, ordinations based on bacterial genera and AMR genes gave similar results (Figure  396 5b). For samples from Bulgaria, Denmark, Italy and especially Germany however, the Procrustes residuals 397 were larger. This was less evident for poultry, though a single Danish poultry herd had a very unusual 398 resistome, considering its taxonomic composition (Figure 5d). Stress plots for NMDS can be found in 399 Supplementary Figure 9. 400

AMR and drug use association 402
We found that the total country-level veterinary AMU was positively associated with AMR in both pigs and 403 poultry. The AMR abundance increased by 1736-3507 (95% CI, β=2621) FPKM in pigs when the AMU 404 increased by 1 loge unit (36.8% increase in AMU) (Figure 6a) and to a lesser degree in poultry, where the 405 AMR abundance increase by 68 -1330 FPKM (95% CI, β =700) when the AMU increase 1 loge unit (Figure  406 6b). For pigs, the variance between farms within countries was 7 times larger than the variance between 407 countries in general, whereas in poultry the variance was 4 times larger within-country than between 408 countries. 409

416
To test if the AMU pattern was associated with AMR gene profiles, we compared the AMR abundance 417 matrix with AMU matrix, comprised of the 15 recorded AM classes (Supplementary Table 5). We found an 418 association between the veterinary AMU pattern and the pig resistomes (Procrustes correlation: 0.51, 419 p<0.001) (Figure 6cd). The poultry resistomes were also significantly associated with AMU, albeit with a 420 lower correlation, likely due to their larger beta-diversity and lower degree of country clustering 421 (Procrustes correlation: 0.38, p<0.001) (Figure 6ef). 422 423

Functional AMR genes 425
In addition to using ResFinder, we also ran most analyses with the FRD database, to elucidate whether the 426 functionally determined AMR genes behave similarly to the acquired AMR genes in ResFinder. If FRD genes 427 serve similar AMR functionality as the acquired ResFinder genes, we would expect similar results. 428 country spread of total AMR, but perhaps the more heterogeneous production system and production 516 management is responsible. 517 DNA extractions from the pooled poultry samples resulted in relatively low DNA yields. The protocol used 518 was optimized for pig feces, human feces and sewage, but not poultry feces. 18 The lower yields 519 necessitated the use of a PCR-based library preparation kit, that can influence downstream analysis of 520 shotgun sequencing. 42 . While the large difference between pig and poultry resistomes in our study is likely 521 real, we caution the use of sensitive, quantitative analyses when comparing between samples prepared 522 using different library preparation kits. For this reason, we have mostly tested within each reservoir. For 523 future studies, one might consider using larger volumes of bird feces or otherwise optimize the protocol to 524 ensure PCR-free library preparation is possible for all samples. 525 The sensitivity of metagenomic approaches does not yet rival phenotypic alternatives such as selective 526 enrichment. We relied on read-mapping as opposed to a metagenomic assembly-based strategy for greater 527 sensitivity. Still, there are AMR genes in important pathogens that we know are likely present but are below 528 our detection limit. For example, we only found evidence for blaCTX-M in three pig herds, whereas in 529 phenotypic studies, the prevalence is high even among farms with no cephalosporin usage. 43 Another 530 tradeoff with the use of reference-based read mapping is our inability to identify mutated housekeeping 531 genes granting AMR or distinguish between close homologs of the same AMR gene, e.g., the non-ESBL 532 blaTEM genes from those encoding ESBL variants. 533 The primary concern with read-mapping techniques, the lack of genomic context, can be solved using 534 metagenomic assembly and binning approaches. 16,44,45 In this way, AMR alleles in full length, their genomic 535 context and their associated taxa have been identified in both pig, poultry and human fecal samples. 46 As 536 shown previously, the association between AMR and AMU is similar for metagenomics and traditional 537 phenotypic methods, but several aspects make metagenomics an intriguing monitoring tool. 17 The fact that 538 both types of analyses (quantitative, sensitive read mapping and qualitative, context-giving binning) use the 539 same raw data, makes metagenomics an attractive tool. In addition, the digital nature of shotgun data 540 would also allow future re-use and form the basis of an invaluable historical archive, potentially usable for 541 both AMR and pathogen tracking worldwide. While it is possible to identify thousands of AMR genes from 542 environmental samples using functional metagenomics, and then track them using shotgun sequencing, 543 their association to drug use and potential risk to human and animal health requires much further work to 544 inform effective drug policies. 545 We found that the metagenomic resistome varied significantly in livestock, with large differences between 546 the pig and poultry reservoir, but also within each livestock species, in a country-dependent manner. 547 Within each country, we found different levels of variation, with some countries having more homogenous 548 herds than others. Differences were seen both in total AMR abundance, but also abundances of AMR types 549 and specific genes, including clinically relevant AMR genes. Some of this variation, we attributed to 550 differential drug usage between the countries. We also identified the microbiome background as an 551 important factor in determining the resistome in livestock, but found the strength of the association was 552 country-dependent, at least in pigs. Interestingly, we found that AMR richness in one livestock species in a 553 country is linked to the abundance in another livestock species. Finally, we observed some indications that 554 newly described AMR genes from functionally metagenomic studies, might not provide AMR functionality 555 when expressed in their natural host, even though they have the potential at the right expression levels in 556 the right organism. 557