Description and determinants of the faecal resistome and microbiome of farmers and slaughterhouse workers_ A metagenome-wide cross-sectional study

Background: By studying the entire human faecal resistome and associated microbiome, the diversity and abundance of faecal antimicrobial resistance genes (ARGs) can be comprehensively characterized. Prior culturebased studies have shown associations between occupational exposure to livestock and carriage of specific antimicrobial resistant bacteria. Using shotgun metagenomics, the present study investigated 194 faecal resistomes and bacteriomes from humans occupationally exposed to ARGs in livestock (i.e. pig and poultry farmers, employees and family members and pig slaughterhouse workers) and a control population (Lifelines cohort) in the Netherlands. In addition, we sought to identify determinants for the human resistome and bacteriome composition by applying a combination of multivariate (NMDS, PERMANOVA, SIMPER and DESeq2 analysis) and multivariable regression analysis techniques. Results: Pig slaughterhouse workers and pig farmers carried higher total ARG abundances in their stools compared to broiler farmers and control subjects. Tetracycline, β-lactam and macrolide resistance gene clusters dominated the resistome of all studied groups. No significant resistome alpha diversity differences were found among the four populations. However, the resistome beta diversity showed a separation of the mean resistome composition of pig and pork exposed workers from broiler farmers and controls, independent of their antimicrobial use. We demonstrated differences in resistome composition between slaughter line positions, pig https://doi.org/10.1016/j.envint.2020.105939 Received 30 April 2020; Received in revised form 29 June 2020; Accepted 30 June 2020 Abbreviations: 90%ID, Clustering of antimicrobial resistance genes was performed at the 90% identity level; AMR, Antimicrobial Resistance; ANOVA, Analysis of Variance; ARG(s), Antimicrobial Resistance Gene(s); BF(s), Broiler farmer(s), including employees and family members; CI, Confidence Interval; CO(s), Control subject(s); DESeq2, Differential gene expression analysis based on the negative binomial distribution; EFFORT, Ecology from Farm to Fork Of microbial Resistance and Transmission project; ESBL-PE, Extended-Spectrum β-Lactamase Producing Enterobacteriaceae; FDR, False Discovery Rate, Benjamini-Hochberg procedure; FPKM, Fragments Per Kilobase reference per Million bacterial fragments; MRSA, Methicillin-resistant Staphylococcus aureus; NMDS, Non-Metric Multidimensional Scaling; PERMANOVA, Permutational Analysis of Variance; PF(s), Pig farmer(s), including employees and family members; SIMPER, Similarity Percentages analysis; Welch’s test, Welch’s Heteroscedastic F Test With Trimmed Means And Winsorized Variances; WO(s), Pig slaughterhouse worker(s); WT, Pairwise Wilcoxon rank-sum Test ⁎ Corresponding author at: Yalelaan 2, 3584CM Utrecht, the Netherlands. E-mail address: l.vangompel@uu.nl (L. Van Gompel). 1 Current affiliation: Reference and Research Laboratory on Antimicrobial Resistance and Healthcare-associated infections, National Microbiology Center, Institute of Health Carlos III, Carretera de Majadahonda Pozuelo, Km. 2.200. 28220 Majadahonda, Madrid, Spain. Environment International 143 (2020) 105939 0160-4120/ © 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/). T versus poultry exposed workers, as well as differences between farmers and employees versus family members. In addition, we found a significant correlation between the bacteriome and resistome, and significant differences in the bacteriome composition between and within the studied subpopulations. Finally, an in-depth analysis of pig and poultry farms – of which also farm livestock resistomes were analysed – showed positive associations between the number of on-farm working hours and human faecal AMR loads. Conclusion: We found that the total normalized faecal ARG carriage was larger in persons working in the Dutch pork production chain compared to poultry farmers and controls. Additionally, we showed significant differences in resistome and bacteriome composition of pig and pork exposed workers compared to a control group, as well as within-population (farms, slaughterhouse) compositional differences. The number of on-farm working hours and the farm type (pig or broiler) that persons live or work on are determinants for the human faecal resistome. Overall, our results may suggest direct or indirect livestock contact as a determinant for human ARG carriage. Future studies should further focus on the connection between the human and livestock resistome (i.e. transmission routes) to substantiate the evidence for livestock-associated resistome acquisition.

versus poultry exposed workers, as well as differences between farmers and employees versus family members. In addition, we found a significant correlation between the bacteriome and resistome, and significant differences in the bacteriome composition between and within the studied subpopulations. Finally, an in-depth analysis of pig and poultry farms -of which also farm livestock resistomes were analysed -showed positive associations between the number of on-farm working hours and human faecal AMR loads. Conclusion: We found that the total normalized faecal ARG carriage was larger in persons working in the Dutch pork production chain compared to poultry farmers and controls. Additionally, we showed significant differences in resistome and bacteriome composition of pig and pork exposed workers compared to a control group, as well as within-population (farms, slaughterhouse) compositional differences. The number of on-farm working hours and the farm type (pig or broiler) that persons live or work on are determinants for the human faecal resistome. Overall, our results may suggest direct or indirect livestock contact as a determinant for human ARG carriage. Future studies should further focus on the connection between the human and livestock resistome (i.e. transmission routes) to substantiate the evidence for livestock-associated resistome acquisition.

Introduction
Recent technological advances have resulted in a rapidly increasing number of studies using metagenomic sequencing to investigate the composition and determinants of the human gut resistome (i.e. the full collection of antimicrobial resistance genes (ARGs) in gut microbiota) (D'Costa et al., 2006). Studies investigating determinants for the human faecal resistome have reported an increase in ARG abundances after antimicrobial use (AMU) or international travel (Bengtsson-Palme et al., 2015;Buelow et al., 2014;Pérez-Cobas et al., 2013), whereas large resistome differences have been demonstrated between countries, including significant associations with country-dependent use of human and veterinary antimicrobials (Forslund et al., 2013(Forslund et al., , 2014Hu et al., 2013). Other studies have described associations between the bacterial community composition and the human/animal faecal or soil resistome (Forsberg et al., 2014;Munk et al., 2018;Pehrsson et al., 2016).
Numerous studies have highlighted occupational exposure of humans to AMR bacteria among agricultural or food industry workers (e.g. farmers, slaughterhouse staff or veterinarians) through direct or indirect contact with livestock (Bisdorff et al., 2012;Dohmen et al., 2017aDohmen et al., , 2017bDorado-García et al., 2015;Gilbert et al., 2012;Graveland et al., 2011;Mulders et al., 2010;van Cleef et al., 2014;Wulf et al., 2008;Bos et al., 2016). However, research has often focused on the transmission of specific resistant bacteria like livestock-associated methicillin-resistant Staphylococcus aureus (LA-MRSA) or extended-spectrum β-lactamase (ESBL) producing bacteria. Studying the entire resistome of occupationally exposed populations could provide a more comprehensive characterization of ARG exposure, offering additional insight into the biological pathways of AMR in these groups.
Our recent studies as part of the European project "Ecology from Farm to Fork Of microbial drug Resistance and Transmission" (EFFORT: http:// www.effort-against-amr.eu/), have shown clear compositional faecal resistome differences between pigs and broilers, and were able to report on determinants for these livestock resistomes Munk et al., 2018;Van Gompel et al., 2019). The current study aims to characterize the abundance and diversity of the faecal resistome and microbiome of 194 persons occupationally exposed to AMR livestock (pig slaughterhouse workers, pig and broiler farmers, their employees and family members working or living in the Netherlands) and 'nonexposed control' subjects from the 'Lifelines' cohort. Additionally, we sought to identify (dis)similarities between, and determinants for, the faecal resistome and microbiome composition of these four populations. Finally, in an in-depth analysis on pig and poultry farms, associations were investigated between the faecal resistome of persons living or working on pig and poultry farms and the faecal resistome of livestock present on these farms.

Study population and data collection
Between October 2014 and December 2015, we visited 31 Dutch conventional, non-mixed, farrow-to-finish pig (N = 19) and broiler farms (N = 12) as part of the EFFORT project. Farm selection criteria were as previously described (Munk et al., 2018). Farmers, their family members and employees living (≥1 day/week) or working (≥0.5 h/ week) on these farms, were invited to hand in a stool sample. We selected all stool samples from adults (≥18 years) for DNA extraction and sequencing from the full pool of participants (N = 94) (i) to omit potential age-related faecal microbiome/resistome differences between children and adults and (ii) to make the age distribution of the farming population more comparable to the population of slaughterhouse workers (below). One broiler farmer's stool sample was removed later in the process for technical reasons resulting in a total of 78 samples (pig farms: N = 54, broiler farms: N = 24). Furthermore, in two oneweek sampling rounds (June 2015 and July 2016), 483 human faecal samples were collected at the largest Dutch pig slaughterhouse as described before (Van Gompel et al., 2020). Overall, pig slaughter was divided into two main phases, slaughter in the black area [successively: lairage, stunning, bleeding, scalding/dehairing and carcass scorching (singeing)] and the clean area [successively: evisceration, pluck removal (removal liver, lungs, heart, oesophagus and tongue), carcass splitting, meat inspection, dressings and carcass cooling], followed by a third phase of cutting and deboning. In total, we selected 70 samples for metagenomic sequencing. Sample selection was based on two criteria (i) inclusion of workers positioned at the most important slaughter steps with respect to the expected occupational AMR exposure risk along the slaughter line (high versus low risk) based on the prevailing (AMR/ bacterial meat contamination) risk factor literature while (ii) preserving enough samples per selected slaughter step to maintain statistical contrast (sampled steps are illustrated in Fig. 1C).
All participants (farmers and slaughterhouse staff) filled out a questionnaire regarding personal information and past occupational exposure (e.g. age, gender, smoking, antimicrobial use and animal contact). Antimicrobial use was defined as using any type of antimicrobials in the past year before sampling. When participants had indicated medicine use without any further specification, participants were regarded as potential antimicrobial users. In sensitivity analyses, (potential) antimicrobial users were excluded.
In total, 46 human faecal samples, collected between February 2014 and January 2015 were retrieved through the Dutch Lifelines biobank (www.lifelines.nl) and serve as controls. Lifelines participants were recruited through general practitioners in the northern part of the Netherlands (Groningen, Friesland and Drenthe) (Stolk et al., 2008). Based on the available questionnaire data, we excluded intravenous and oral antimicrobial users (≤three months before sampling) and all individuals who reported to live on a farm, worked with farm animals or reported to be butchers. All Lifelines participants were not hospitalized three months prior to the assessment. All but one participant were born in the Netherlands.
Overall, the farming population consisted of farmers, employees, spouses/life partners, and other family members. Hereafter, we will refer to the overall subpopulations as 'pig farmers, employees and family members' (PFs), 'broiler farmers, employees and family members' (BFs), 'pig slaughterhouse workers' (WOs) and 'controls' (COs). To refer to results from the farming subpopulations, the terms 'farmers', 'employees' and 'family members' ('spouse' and 'other family') will specifically be used.

Sample handling, DNA extraction and first laboratory preparation steps
Farmers were instructed to collect their faecal sample as close as possible to the collection by the researchers and to keep the sample frozen. Following the farm visit, samples were transported on dry ice and stored at −80°C the same day.  Fig. 1. Overview of the normalized faecal AMR abundances per and within study subpopulation(s). 1a: Overall ARG abundance (FPKM) per study subpopulation. 1b: Overall ARG abundance (FPKM) in the farming subpopulations. 1c: Overall ARG abundance (FPKM) per slaughter position in the slaughterhouse population. Boxplots: A box represents the 25th (Q1) to 75th (Q3) percentile, the centreline depicts the median (Q2). Dots represent values larger or smaller than Q1 -1.5 or Q3 + 1.5 times the interquartile range. N = Number of observations in each group. *1a: 'pig sh worker' = pig slaughterhouse worker. 'pig farmer' and 'broiler farmer' categories include employees and family members. The Welch's Heteroscedastic F Test with trimmed means and Winsorized Variances with trimming rate 0.1 tests differences between the four subpopulations. Values above the arches represent the p-values corresponding to the pairwise tests adjusted for multiple testing (Benjamini-Hochberg procedure). *1b: The Kruskal-Wallis test tests farmers and employees versus family members (the latter include spouses and other family members). *1c: The Kruskal-Wallis test includes all slaughter positions presented on the X-axis. The slaughter line runs from left to right. The cooling area is located between the job position 'other removal' and the 'deboning' room.
were collected and stored as described before (Van Gompel et al., 2020). Lifelines' faecal samples were kept refrigerated by the participant until collection by the researchers, followed by transportation on dry ice and storage at -80°C until DNA extraction.
Samples were thawed on ice and of each sample, DNA was extracted using the modified QIAamp® Fast DNA stool mini kit (Qiagen, cat. no. 51604) as previously applied and demonstrated to preserve microbiome diversity (Knudsen et al., 2016;Munk et al., 2018). DNA was subsequently stored at −80°C until further processing. Frozen DNA was consecutively shipped on dry ice to Wageningen Bioveterinary Research (Lelystad, the Netherlands) for library preparation. Libraries were prepared with the amplification poor NEBNext® Ultra DNA Library preparation kit (3-10 cycles, New England BioLabs®Inc., UK). Multiplexing and sequencing were performed using the Illumina HiSeq4000 platform (2x150bp paired-end sequencing) targeted at 40 million PE clusters (80 M reads) per sample (GenomeScan, Leiden, the Netherlands).

Metagenomic processing and analysis
Bioinformatic processing was performed as described before (Munk et al., 2018). Briefly, quality controlled and adapter-cleaned read-pairs were mapped to the NCBI bacterial genome reference database and the ResFinder reference resistome database. Raw counts matrices were preserved for certain analyses (see below), while the bacteriome count matrix was corrected for genome length. To normalize the resistome counts, Fragments Per Kilobase reference per Million bacterial fragments (FPKM) were computed. The bacteriome was classified at the taxonomic phylum, class and genus level, whereas the resistome was clustered at the AMR class and the 90% identity cluster level (ARG cluster, CD-HIT-EST) (Fu et al., 2012;Munk et al., 2018).

Database management and statistical analysis
Questionnaires were entered and checked for consistency in EpiData 3.1 before importing into 'R' for data cleaning and statistical analysis (an overview of the used 'R' version and packages, including references for each package, is provided in the supplement). The most important potential ARG determinants were selected from the questionnaire for further use in the analysis when at least 95% of the questions were answered and if the question still provided enough contrast. Slaughterhouse participants were asked to indicate the slaughterhouse position(s) where they worked the most (maximum of three) from a list of slaughter line positions. In case workers had indicated multiple job positions, slaughterhouse staff was assigned to the last position in the slaughter line as previously described (Van Gompel et al., 2020). For all analyses, we initially grouped and compared the four study populations based on the expected intensity of animal contact (PFs > BFs > COs), animal species (pigs versus poultry), and the occupational setting (farms versus slaughterhouse), followed by within-population analyses (farms, slaughterhouse).

Alpha diversity
Overall richness (observed number of species) and the Shannon diversity index (a combination of richness and evenness), both alpha diversity indices, were calculated after rarefying the raw resistome count matrix and the genome length corrected species matrix (multiplied by 10 6 before transforming to integers) while retaining at least 95% of all samples (90% per group) ('R' vegan function 'rrarefy'). Subsequently, overall statistical differences and pairwise comparisons between the different populations were computed followed by adjustments for multiple testing [False Discovery Rate (FDR), Benjamini-Hochberg procedure].

Clustering of and identifying determinants for the human faecal resistome and bacteriome (NMDS, PERMANOVA, Procrustes analysis)
Before ordinating the faecal resistome (FPKM, at AMR class and ARG cluster level) and bacteriome (genome length corrected and normalized by the sum of genome length corrected counts at the genus level), Bray-Curtis (BC) dissimilarities were calculated. Non-metric multidimensional scaling (NMDS), an indirect gradient analysis, was performed in a 2-dimensional configuration ('R' vegan function 'metaMDS'). In case of stress levels > 0.15, NMDS was repeated in three dimensions. Additionally, differences in overall resistome or bacteriome composition and/or relative abundances between groups were tested within a permutational multivariate analysis of variance (PERMAN-OVA, BC-dissimilarities, 'R' vegan function 'adonis' at 999 permutations) (Anderson, 2001). PERMANOVA was applied to identify possible determinants for between and within-population resistome and bacteriome dissimilarities. For each PERMANOVA test, the assumption of homogeneity of variances was checked (p > 0.05) using the 'R' vegan function 'betadisper'.
The association between the resistome (90% ID clustering) and bacteriome (genus level) was determined by superimposing both resistome and bacteriome NMDS outcomes (2D) in a symmetric Procrustes analysis ('R' vegan function 'procrustes'), an analysis in which both ordinations ('shapes') are compared to one another retrieving a goodness-of-fit statistic (Munk et al., 2018). A p-value was retrieved through the 'R' vegan function 'protest' at 999 permutations.

Bray-Curtis dissimilarity decomposition (SIMPER) and differential abundance analysis (DESeq2)
The most interesting resistome or bacteriome pattern of dissimilarity identified within the NMDS and PERMANOVA analyses was further investigated by means of a similarity percentages analysis (SIMPER, 'R' vegan function 'simper') to assess the contribution of the individual ARGs or bacterial genera to the overall dissimilarities displayed in the beta diversity ordinations (Clarke, 1993). SIMPER uses the input as described for the NMDS, and defines the contribution of each ARG cluster in explaining at least 70% of the overall dissimilarities between populations as visualized within the NMDS (resulting from variation in feature abundances and differences among groups). Specific differences between groups were further detailed within a differential abundance analysis (Bioconductor differential gene expression analysis based on the negative binomial distribution: DESeq2 'R' package). DESeq2 uses empirical Bayesian shrinkage for dispersion estimation and a negative binomial distribution to estimate log fold changes of features between groups (Love et al., 2014). For our study, we used the raw resistome data, and applied no further shrinkage for log fold change estimation, while size factors (total bacterial pairs divided by the minimum number of total bacterial pairs) (Munk et al., 2018) were provided for normalization purposes. FDR corrected p-values were obtained as previously described (Munk et al., 2018).

In-depth analysis of the faecal resistomes of pig and broiler farmers and their livestock
The faecal resistomes of the PFs and BFs were analysed together with the previously reported faecal resistomes of pigs and broilers sampled at the corresponding Dutch farms. Per farm, 25 faecal samples from pigs and poultry were previously collected, pooled, sequenced and metagenomically analyzed as described before (Munk et al., 2018). Resistome counts were expressed as FPKM values. DNA sequences of these samples were deposited at the European Nucleotide Archive (project accession number: PRJEB22062).
In total, five multivariable regression models were fitted: each of the four highest abundant human faecal AMR class clusters (i.e. tetracycline, β-lactam, macrolide, aminoglycoside) and the overall AMR load detected in human stool, were separately regressed against their respective AMR class clusters found in pig and poultry faeces. All models were adjusted for potential confounding: i.e. gender, personal antimicrobial use, current smoking, age, number of on-farm working hours, type of person (farmer and employee versus family members) and type of farm (living/working on a pig versus poultry farm). First, the full models were fitted and compared with or without a random intercept for farm (linear mixed models: 'R' lme4 function 'lmer', linear regression: stats function 'lm'). Variables were log e transformed if deemed necessary. Then, the models were reduced manually (mixed models: Satterthwaite's method of single term deletions and the maximum likelihood method) or automatically (linear regression: stepwise backward/forward selection, 'R' stats function 'step'), after which the model assumptions were checked. Due to the small sample size, the models were fitted for both the separate and the combined pig and broiler 'human-animal' datasets.
A possible correlation between the human resistome and the animal resistome was assessed using a Procrustes analysis (procedure as described in Section 2.6., NMDS in two dimensions). To match the number of livestock samples in the Procrustes analysis, the human resistome was downsampled to one person per farm by (i) including only the farmer (or if no farmer was present the employee, and if also no employee was present a family member) or (ii) the person with the highest number of on-farm working hours (if > 1 person remained per farm, the farmer was selected).
Finally, overall human-livestock resistome compositions were visualized in a heatmap. The heatmap includes all ARGs which are present in at least one livestock and one human sample. Abundances were log-transformed [log e (ARG + 1)]. One dendrogram compares the human and animal samples based on the Bray-Curtis dissimilarity matrix (untransformed ARG abundances, hierarchical clustering). ARGs (second dendrogram) are clustered based on the Spearman correlation coefficients.
Participants had an average age of 42.5 years and were mostly nonsmokers (71.7%). Within the occupationally exposed populations (farming populations, pig slaughterhouse workers), between 8.6 and 13% had used antimicrobials of any kind in the past year before sampling. Within the farming subpopulation, the largest groups were farmers (38.5%) and their spouses (33.3%). All farmers and employees entered the stables and nearly all had daily contact with their farm animals (≥ 92.6%). While the majority of the pig and broiler farmer family members did enter the stables (≥ 83.3%), the majority had no daily animal contact (≥ 66.7%) (Supplementary Table S1). The mean reported number of working hours per week on pig farms (32.6 h) was higher than on broiler farms (23.5 h), but only when family members were excluded a (borderline) significant difference (Analysis of Variance (ANOVA), p = 0.05) was noted between the mean number of hours worked by pig farmers and employees (48.5 h) compared to broiler farmers and employees (37.7 h) (Supplementary Table S1).
Faecal resistome counts were normalized for bacterial counts (FPKM values) and clustered at two levels [90% identity level (ARGs) and the AMR phenotypic drug class level (AMR class cluster level) as defined in the ResFinder database]. Absolute levels of ARG abundance (total FPKM) were significantly higher in both PFs and WOs compared to BFs and COs (Fig. 1A). At the AMR class level, we observed large median differences between these populations for e.g. β-lactam and macrolide resistance clusters (Supplementary Fig. S1A). At pig farms, the highest median ARG abundances were seen in farmers and employees (Fig. 1B). At broiler farms, the highest median ARG abundances were found in other family members, closely followed by the farmers. When examining persons from both farming subpopulations together, we found a significant difference in median ARG abundance between farmers and employees compared to family members (Kruskal-Wallis test, p = 0.008).
In WOs, a decreasing mean ARG trend is visible along the slaughter line, whereas the highest mean ARG abundances were found at lairage, evisceration and pluck removal. Significantly different total ARG abundances were found between the slaughter line positions presented in Fig. 1C (Kruskal-Wallis, p = 0.04). Additionally, total ARG loads varied significantly between the black area (positions lairage to scalding & dehairing) and the deboning area [Pairwise Wilcoxon ranksum test (WT), p = 0.006], and more specifically between lairage and deboning (WT, p = 0.048). At the AMR class level, a decreasing trend was mainly observed within the (most abundant) tetracycline and βlactam resistance clusters ( Supplementary Fig. S2A-B). Furthermore, a more fluctuating (up and down) trend along the slaughter line was identified for e.g. the macrolide and aminoglycoside resistance clusters. Relative resistome and bacteriome abundances across the slaughter line are further illustrated in the supplement ( Supplementary Fig. S3A-D). Missings are excluded from the summary statistics. N = Number of participants in each group for which data concerning age was present. *One observation was missing. **No antimicrobial use: Participants reported no antimicrobial use within one year before sampling. ***No antimicrobial use: Participants reported no oral or intravenous antimicrobial use up to three months before sample collection.

Faecal resistome and bacteriome richness and evenness
To compute alpha diversity indices, count matrices were rarefied at respectively 5,185 (resistome) and 104,737,719 (bacteriome) counts per sample, retaining respectively 97.4% and 98.5% of all samples. Rarefaction curves are provided in the supplement (Supplementary Fig.  S4). In total, across all human subpopulations, we identified 208 ARG clusters (at 90%ID) of 14 different AMR classes (13 AMR classes in the occupationally exposed populations) of which 38.5% were shared among all subpopulations ( Supplementary Fig. S5). No significant differences in resistome alpha diversity were found between the four study subpopulations (Fig. 3). In contrast, significant overall differences in mean bacteriome alpha diversity were demonstrated (ANOVA/Welch's test: Richness and Shannon diversity index, p ≤ 0.0007). Species richness and Shannon diversity of WO's differed significantly from the three other study subpopulations, while the bacteriome Shannon diversity index of the PFs differed significantly from the COs.

Faecal resistome and bacteriome clustering
Clusters were visualized using Non-Metric Multidimensional Scaling (NMDS). Fig. 4 shows NMDS configurations based on BC dissimilarities for the faecal resistome and bacteriome. The distribution of the resistome population mean (centroid) associated with the BFs and the COs was fairly distinct from those associated with PFs and WOs. Clustering patterns did not largely change after excluding potential antimicrobial users, nor after the exclusion of family members or deboning staff (potential 'low-risk groups') from respectively the farming and slaughterhouse subpopulations. The latter resulted in a slight overlap of the PF and BF's 95% centroid confidence interval (CI, data not shown). Similar compositional differences were observed for the bacteriome, although the population mean (95% CI) of persons sampled at broiler farms was less distinguishable from persons sampled at pig farms (Fig. 4). Additionally excluding above mentioned lower risk categories resulted in overlapping centroids of pig and pork exposed persons (farmers, slaughterhouse staff) (data not shown).

Fig. 2.
Relative faecal resistome and bacteriome abundance per study subpopulation. 2ab: Relative resistome abundances (FPKM) are clustered at AMR class (a) and at the ARG cluster level (90%ID) (b) and are scaled to 100%. 2 cd: Bacteriome phylum (c) and genus (d) level abundances are corrected for genome length and scaled to 100%. 2bcd: Plots represent the top ten most abundant features across the four study subpopulations and labels are arranged in decreasing order with respect to overall ARG abundance across the overall population. All remaining features are part of the category 'Other'. Legend labels are arranged with respect to a decreasing total ARG abundance calculated for the overall study population. X-axis: 'pig sh worker' = pig slaughterhouse worker (WO); 'pig farmer' (PF) and 'broiler farmer' (BF) categories include employees and family members. ARG cluster names: oxaz-phen = oxazolidinone & amphenicol, macr-oxaz-phen = macrolide & oxazolidinone & amphenicol.

Determinants for the faecal resistome and bacteriome of humans occupationally exposed to livestock
Compositional differences between the studied subpopulations Permutational Analysis Of Variance (PERMANOVA) revealed significant differences between the faecal resistome compositions of the four study subpopulations (ARG cluster: R 2 = 12.3%, Table 2; AMR class: R 2 = 16.0%, data not shown, both p = 0.001).
We also identified significant within-population dissimilarities ( Table 2). In particular, small resistome differences were demonstrated distinguishing slaughterhouse staff based on job position [working at the (i) black area, (ii) clean area and (iii) deboning room], which were not found for the bacteriome. Within the farming subpopulation, resistome and bacteriome compositional differences were observed between PFs and BFs, and between farmers & employees and family members.

The human faecal bacteriome shapes the resistome
The Procrustes analysis revealed a significant correlation (0.49, p = 0.001) between the faecal resistome and bacteriome recovered from the overall study population (Supplementary Fig. S6A). The largest range in (dis)similarities between the resistome and bacteriome was seen in the slaughterhouse population ( Supplementary Fig. S6B).

Occupational associations with specific ARGs
NMDS ordinations mainly showed a shift of the mean resistome and bacteriome composition in subjects occupationally exposed to pigs and pork compared to the COs. Therefore, we further examined the compositional dissimilarities between these populations [pig and pork exposed workers (PFs/WOs) versus the control subjects (COs)] using two types of analyses: Similarity Percentages analysis (SIMPER) and Fig. 3. Alpha diversity indices per study subpopulation (ARG richness and Shannon diversity index). Alpha diversity indices were calculated after rarefying the raw resistome count and genome length corrected bacteriome matrices. 3a: Resistome richness. 3b: Resistome Shannon diversity index. 3c: Bacteriome richness. 3d: Bacteriome Shannon diversity index. 'pig sh worker' = pig slaughterhouse worker, 'pig farmer' and 'broiler farmer' categories include employees and family members. Values above the arches represent the p-values corresponding to the pairwise tests adjusted for multiple testing (Benjamini-Hochberg procedure). Welch's test = Welch's Heteroscedastic F Test with trimmed means and Winsorized Variances with trimming rate 0.1. NA = Pairwise comparisons are not applicable. N = Number of observations in each group after rarefaction. Boxplots: A box represents the 25th (Q1) to 75th (Q3) percentile, the centreline depicts the median (Q2). Dots represent values larger or smaller than Q1 -1.5 or Q3 + 1.5 times the interquartile range. DESeq2.

Results Similarity Percentages analysis (SIMPER)
Highly abundant AMR class clusters, tetracycline [tet(Q), tet(W)] and β-lactam (cfxA6, cfxA_clust) resistance clusters, contributed most to the overall resistome dissimilarity between PFs/WOs and the COs as visualized in the NMDS (up to a~64% cumulative contribution, Supplementary Excel Table 7a). The overall bacteriome dissimilarity is mainly resulting from the joint contributions of the following bacterial genera in decreasing order (up to a~70% cumulative contribution): Prevotella, Bacteroides, Cellulomonas, Oscillobacter, Roseburia, Blautia and Alistipes (Supplementary Excel Table 7b).

Results differential abundance analysis (DESeq2)
A more in-depth ARG differential abundance analysis between the pig and pork exposed workers and the control population was performed using DESeq2. Here, we found that resistomes of PFs/WOs contained significantly higher levels of 30 ARGs compared to the resistomes of COs (range: 1.9-32.2 times as high, FDR p < 0.05), with only one ARG cluster (VanR-D) being significantly less abundant in pig/ pork workers compared to controls (Supplementary Excel Table 8). Of these ARGs, 22 genes (73.3%) are part of three highly abundant AMR classes previously identified (tetracycline, β-lactam, macrolide resistance).

Determinants of the resistome of persons working or living on pig and broiler farms
Matching livestock and human faecal samples were available for 76 human samples. A brief description of the resistome and bacteriome of the human and livestock populations is provided in the supplement, as well as a heatmap showing clear separate human and livestock resistome clusters (Supplementary Figs. 8 and 9). After removing variables based on incomplete questionnaire data, a modelling dataset of 73 samples remained which originated from 18 pig farms and 12 broiler farms. Results from the regression analyses including the most common AMR class clusters show that livestock ARG abundances do not have a significant effect on ARG abundance in humans. However, the number of on-farm working hours was positively associated with certain types of AMR found in human stool (i.e. total faecal AMR load, tetracycline, macrolide and aminoglycoside resistance abundance, Supplementary Tables 2-6). Additionally, people living or working on pig farms have significantly higher AMR loads in their stool than persons living or working on broiler farms (i.e. total AMR abundance and β-lactam resistance in the reduced models). Furthermore, farmers and employees tend to carry higher β-lactam, but lower aminoglycoside resistance gene concentrations in their stools compared to family members of farmers. These associations could not always be confirmed when tested in the two farming subpopulations separately, likely in part due to the lower sample size, in particular of the broiler dataset.
Non-significant correlations were found in the Procrustes analysis between the resistome of farmers, workers and family members and the resistome of their livestock when analyzing the pig and broiler farm datasets separately. On broiler farms, the correlation was slightly higher when per farm the principal farmer was selected for the analysis, compared to the correlation in which the person with the highest number of on-farm working hours was selected, whereas in pig farms the correlation more or less remained the same [pig farms: 0.34 (working hours), 0.35 (principal farmer); broiler farms: 0.38 (working hours), 0.33 (principal farmer); all p > 0.1]. When the pig and broiler datasets were analysed together, relatively higher and significant correlations were found [0.48, p = 0.003 (working hours)/0.39, p = 0.01 (principal farmer)].

Discussion
The present study used metagenomic shotgun sequencing to investigate the abundance and diversity of the faecal resistome and bacteriome of persons occupationally exposed to pigs or poultry in the Netherlands. Our goal was to identify determinants for the composition A B Fig. 4. Resistome and bacteriome ordinations (NMDS, Bray-Curtis dissimilarity). 4a: Resistome clustering at ARG cluster level (90%ID). 4b: Bacteriome clustering at the bacterial genus level. 'pig sh worker' = pig slaughterhouse worker. 'pig farmer' and 'broiler farmer' categories include employees and family members. 2D/ 3D = NMDS in 2 or 3 dimensions. stress = stress associated with the NMDS ordination. Ellipses represent the 95% confidence intervals of the centroid's standard error per study subpopulation. The associated Shepard stress plots can be found in Supplementary Fig. S7.
of these resistomes and bacteriomes. Findings from this study show an increased ARG carriage in pig and pork exposed populations (PFs/WOs) compared to control subjects (COs). No such differences were found between broiler farmers (BFs) and COs.

Animal exposure as a potential determinant for resistome composition in humans
Overall, higher ARG loads were identified in PFs and WOs compared to BFs and COs. This difference between PFs and BFs was also confirmed in the regression analyses for overall AMR and β-lactam loads. Furthermore, we did not observe significant differences in resistome alpha diversity among the four studied subpopulations, but showed differences in mean resistome composition between pig and pork exposed workers and family members (PFs/WOs) and BFs/COs in the ordinations. Using PERMANOVA, we also observed significant resistome compositional differences between these populations. The effect was shown regardless of antimicrobial use, a potential determinant of faecal microbiome composition and ARG carriage (Buelow et al., 2014;Dethlefsen and Relman 2011;Kirchner et al., 2014). In addition, removing potentially lower-exposed persons (farmers' family members, slaughterhouse 'deboning' staff) from the analysis increased the explained variation substantially. Further comparison of the resistome of PFs/WOs with the COs revealed higher mean ARG abundances within faecal samples from PFs/WOs compared to those of COs for a limited number of (mainly) highly abundant ARG clusters. Finally, next to between-subpopulation differences, we also demonstrated smaller significant within-subpopulation resistome compositional differences based on farm type (pig versus broiler), type of person (farmers and employees versus family members), and job position within the slaughterhouse.
To the best of our knowledge, no direct comparisons between persons sampled in slaughterhouses, farms or the general population are currently available regarding resistome composition or overall ARG loads. Prior culture-based studies in occupationally exposed populations have shown increased carriage rates of specific resistant bacteria [nasal LA-MRSA, faecal ESBL-producing Enterobacteriaceae (ESBL-PE)]. For instance, an overall higher MRSA prevalence was reported in pig farmers compared to pig slaughterhouse staff or broiler farmers, with the two latter populations in turn carrying a higher prevalence compared to the general population (Geenen et al., 2013;Gilbert et al., 2012;Mulders et al., 2010;Reynaga et al., 2016;van Cleef et al., 2010). This contrasts with reports on overall faecal ESBL carriage, which was described to be higher in broiler farmers compared to pig farmers and pig slaughterhouse workers employed at early positions in the pig slaughter line, while the latter two populations again carried similar faecal ESBL loads in comparison to the general population (Dohmen et al., , 2017bHuijbers et al., 2014). Possible explanations for these different carrier rates among occupationally exposed populations include differences in bacterial prevalences found between livestock species to which farmers or slaughterhouse workers are exposed to. Resistome compositional differences have been reported between pigs and poultry, while higher mean ARG loads were identified in pigs Table 2 Faecal resistome and bacteriome compositional differences between occupational exposed groups incl. and excl. antimicrobial users. compared to broilers (Munk et al., 2018), levels which were confirmed for the pooled pig and broiler faecal samples from the Dutch farms included in the current study: the mean faecal ARG load in pigs was 4352 FPKM (10-90th percentile: 3084-6472 FPKM) versus a mean ARG load of 1880 FPKM in broilers (10-90th percentile: 1061-2814 FPKM). Culture-based studies have also described positive associations between nasal or faecal AMR carriage and on-farm working hours, intensive on-farm livestock contact (especially with MRSA or ESBL-positive livestock), and working in the stables Dorado-García et al., 2015;Huijbers et al., 2014;van Cleef et al., 2014van Cleef et al., , 2015Van Den Broek et al., 2009). Similarly, having contact with live animals or working at early slaughter positions in the slaughter line have been associated with AMR carriage in slaughterhouse staff (Dohmen et al., 2017b;Gilbert et al., 2012;van Cleef et al., 2010). The regression analyses demonstrated significant associations between the number of on-farm working hours and/or the person sampled (farmers and employees versus family members), and the concentration of certain types of AMR found in human stool. However, no association was found between AMR present in livestock faeces and AMR present in the stool of FPs and BFs, while significant human-livestock resistome correlations were found only when pig and broiler datasets were combined in the analysis (probably caused by the strong clustering per animal species). Therefore, direct animal exposure only is unlikely to explain the differences in ARG abundances or resistome dissimilarities between the pig and broiler farm subpopulations. It could be hypothesized that indirect environmental exposure to AMR bacteria (e.g. transmission through air) might result in AMR carriage. Accordingly, prior culturebased AMR studies have shown associations between AMR bacteria in barn dust and human AMR carriage (Dohmen et al., 2017a;Bos et al., 2016). The overall potential for direct or indirect bacterial and/or ARG transmission might also be higher at pig farms than broiler farms. Although we did not measure the intensity of individual animal contact in this study (e.g. recording of specific on-farm activities), conventional Dutch farrow-to-finish farmers in the Netherlands tend to engage in a higher degree of individual animal contact compared to conventional Dutch broiler farmers, as they are involved in various activities such as inseminating, handling of piglets (e.g. applying ear tags, tail cutting), weaning piglets or individual medicine treatments, activities which are not part of conventional broiler farming. Since we did find a positive association between ARG carriage and on-farm working hours, but not between human ARG carriage and ARG carriage in livestock, it could also be speculated that not the most optimal gene combinations were tested in our regression models. A link between the livestock and human resistome might e.g. only be found for those resistome markers which permit distinguishing between livestock reservoirs [to be defined in source attribution studies (Ribeiro Duarte et al., data in preparation)].

Resistome differences within the context of the slaughterhouse
We demonstrated relative resistome dissimilarities between pig slaughterhouse staff. In particular, we identified a decreasing trend in overall faecal ARG loads along the slaughter line, with significantly higher overall ARG loads in staff working at earlier slaughter line positions (black area, lairage) compared to persons working at the deboning room. Additionally, we found significant resistome compositional differences with regard to job location in the slaughterhouse [(i) black, (ii) clean, or (iii) deboning area]. These findings confirm results from prior culture-based research in pig slaughterhouse employees, where higher nasal MRSA and faecal ESBL carriage were found to be associated with live animal contact or working at early slaughter line positions (Dohmen et al., 2017b;Gilbert et al., 2012). Our recent qPCR study tested associations between slaughter line job position and tet (W)/erm(B) human environmental exposure to these ARGs and human faecal tet(W)/erm(B) carriage (Van Gompel et al., 2020). The latter study identified clear decreasing ARG loads [tet(W) and erm(B)] in the environment along the slaughter line, but no robust association was found between environmental ARG exposure and human faecal ARG carriage.

The role of the bacteriome
Prior research has highlighted the importance of the microbiome in shaping the resistome (Forsberg et al., 2014;Pehrsson et al., 2016). Correspondingly, we found a significant correlation between bacteriome and resistome compositions of the overall population. Although the current study design does not allow for analysis of the actual resistome/bacteriome dynamics in detail, the observed correlation between the resistome and bacteriome datasets emphasises potential bacteriome/resistome intertwinements and stresses the importance of considering the bacteriome when analysing resistance determinants.
Significant differences between bacteriome alpha diversity of WOs and other subpopulations were demonstrated, as well as a significant difference between bacteriome alpha diversity in PFs and COs (mainly due to evenness: linearly extrapolated from a significant difference in the Shannon index and a non-significant difference in ARG richness). Similar to the resistome, bacteriome compositions of pig and pork exposed populations were shifted from the controls. In addition, we found within-farm bacteriome dissimilarities depending on farm type or type of person (farmers and employees versus family members). Interestingly, the highest explained bacteriome variation was found when excluding potentially lower-exposed family members and slaughterhouse deboning staff, which seems to suggest an effect of livestock exposure.
The healthy microbiome is known to be heavily influenced by environmental factors like lifestyle and diet (David et al., 2014;Rampelli et al., 2015;Rothschild et al., 2018). The majority (~85%) of the slaughterhouse population consisted of Eastern European workers (mostly Polish and Romanian nationalities), often temporarily residing in the Netherlands, which are potentially exposed to different diets, habits and environments (e.g. through long-distance commuting). The latter might explain differences in gene richness compared to the other three farming and control populations assumed to be Dutch or more permanently residing in the Netherlands. To the best of our knowledge, the role of frequent and/or close human-livestock contact in shaping the adult gut microbiome has not been investigated before, yet, exposure to pets was previously shown to influence the infant gut microbiome (Tun et al., 2017). Others have suggested that factors such as cohabitation with dogs, living in urban versus rural areas, the type of occupation (dairy farmers versus non-farmers) and contact with pig barn indoor-air could affect the human skin or nasal microbiome (Kraemer et al., 2019;Shukla et al., 2017;Song et al., 2013;Ying et al., 2015).

Overview of future study improvements
Our cross-sectional study provides first insights into the composition and determinants of the faecal resistome and microbiome of populations occupationally exposed to ARGs in livestock. In contrast to previous culture-based studies, we did not study ARG expression in viable (pathogenic) bacteria, but instead provided a comprehensive overview of ARGs in total community DNA. Within our study, we shotgun sequenced and analysed a relatively large dataset of nearly 200 human faecal samples. The increasingly lower costs of sequencing and the continuous development of improved analysis techniques will provide ample future opportunities to further examine other large AMR exposed populations and additional subpopulations (e.g. within the slaughterhouse), preferably using a longitudinal design or an assembly-based metagenomic strategy to better understand cause-effect and microbiome-resistome relationships.
Based on our study results, we hypothesize that the faecal resistome of persons occupationally exposed to livestock might be influenced by direct contact with livestock or indirectly through contact with the farming or slaughterhouse environment [e.g. indirect contact with ARGs and AMR bacteria on carcasses, in barn dust and air]. To further substantiate the former claim, it is recommended to include larger livestock resistome datasets or other environmental data into the models. Pig and broiler farm dust sampling and sequencing were included in the EFFORT project and are further analysed as part of a follow-up study in which the importance of the farm environment is highlighted (Luiken et al., data in preparation).

Conclusion
We identified higher faecal ARG loads in pig and pork exposed workers compared to broiler farmers and control subjects, as well as significant faecal resistome and bacteriome dissimilarities between and within (farms, slaughterhouse) populations occupationally exposed to AMR. On-farm working hours and working or living on a pig farm (versus broiler farm) are determinants for the resistome of people living or working on farms. Our results suggest that livestock contact may be a determinant for faecal ARG carriage. Further studies are needed to examine the influence of (in)direct livestock contact on the human faecal resistome.

Ethics approval and consent to participate
The Medical Ethical Committee of the University Medical Centre Utrecht (NL) confirmed that the Dutch 'Medical Research Involving Human Subjects Act' did not apply for the study of the EFFORT populations (Protocols 14-346/C, 14-403/C). 'Lifelines' research operated according to the protocols approved by the Medical Ethics Review Board of the Medical Center Groningen (NL) (Protocol METc2007/ 152). All participants gave written consent. EFFORT participants were financially compensated (25 euro: slaughterhouse employees, 100 euro: per farm family which includes efforts for participation in previous studies (Munk et al., 2018)).

Data disclosure and availability of data and materials
This research uses datasets from three different study populations. Participants from the farming population (EFFORT project) were selected from Dutch pig and poultry farms sampled within the scope of Munk et al. (2018). The slaughterhouse population (EFFORT project) was partly analysed for ESBL in Dohmen et al., (Dohmen et al., 2017b) (cohort 2015) and tested for specific ARGs (tetW and ermB, cohort 2015 and 2016) in Van Gompel et al., (Van Gompel et al., 2020). The Dutch control population is selected from the Dutch Lifelines cohort (Stolk et al., 2008). The 194 resistomes and bacteriomes of the farmers, slaughterhouse workers and the Lifelines control subjects, analysed and described in this paper, have not been published before.
The DNA sequences (reads) from the 194 metagenomic samples are deposited in the European Nucleotide Archive (EGA) under project accession number EGAS00001003944. All data access should follow the informed consent regulations as stipulated within EFFORT and Lifelines. Access to the metadata from the control group ('Lifelines' cohort: e.g. age, gender, antimicrobial use, animal contact) was purchased from Lifelines for a period of 6 months (data access agreement OV19_0483) and can only be retrieved through www.lifelines.nl/researcher. This excludes the variable gender which can be retrieved through the EGA repository. Nevertheless, Lifelines metadata was only used to select a suitable control group and was not used in the analyses. Finally, due to the potential of tracing back individuals at the slaughterhouse based on the combination of job position, age and nationality, open-access metadata from this group is restricted to the variables directly used in the main analysis (job position, antimicrobial use) and the additional variable sampling round.

Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Vion Food Group (owner of the slaughterhouse in which workers were sampled) was not involved in the statistical analysis, the interpretation of the data, nor the decision to publish.
Netherlands. It employs a broad range of investigative procedures in assessing the biomedical, socio-demographic, behavioural, physical and psychological factors which contribute to the health and the disease of the general population, with a special focus on multi-morbidity and complex genetics.