Machine learning and metagenomics reveal shared antimicrobial resistance profiles across multiple chicken farms and abattoirs in China

Abstract

China is the largest global consumer of antimicrobials and improving surveillance methods could help to reduce antimicrobial resistance (AMR) spread. Here we report the surveillance of ten large-scale chicken farms and four connected abattoirs in three Chinese provinces over 2.5 years. Using a data mining approach based on machine learning, we analysed 461 microbiomes from birds, carcasses and environments, identifying 145 potentially mobile antibiotic resistance genes (ARGs) shared between chickens and environments across all farms. A core set of 233 ARGs and 186 microbial species extracted from the chicken gut microbiome correlated with the AMR profiles of Escherichia coli colonizing the same gut, including Arcobacter, Acinetobacter and Sphingobacterium, clinically relevant for humans, and 38 clinically relevant ARGs. Temperature and humidity in the barns were also correlated with ARG presence. We reveal an intricate network of correlations between environments, microbial communities and AMR, suggesting multiple routes to improving AMR surveillance in livestock production.

Fig. 1. Analysis of potentially mobile ARGs.

a, Pie chart showing the proportion of ARGs (out of the 195 found) associated with one or multiple MGEs. b, Undirected network graph showing potentially mobile ARGs (small orange circles) associated with different sample sources (large green circles). The edges in the graph link the potentially mobile ARGs to the sources in which they were found. c, Number of potentially mobile ARGs per sample per source. Each circle represents a single sample, with circles coloured by farm. d, Venn diagram showing that 145 (out of 661) potentially mobile ARGs were found to be present in both chicken and environmental samples from the same farm, and 182 potentially mobile ARGs contained clinically relevant ARGs. An overlap of 46 clinically relevant potentially mobile ARGs was found in chicken and environmental sources obtained from the same farm. Note that in this analysis, samples from the same source collected at t₁ (week 3) and t₂ (week 6) were aggregated together, leading to a total of seven sources considered for each farm.

Fig. 2. Data mining pipeline to find correlations between gut microbiome, antibiotic resistance in E. coli, temperature and humidity.

The full data analysis workflow of the bespoke data mining method based on ML. Input data are shown in green. Phase I involves metagenome data pre-processing (in yellow). The steps are described in detail in the Methods section. Phase II involves the training and testing of ML-powered predictive functions to isolate metagenomic features (that is, the ARG count and relative abundances of microbial species present in the sample) correlated with phenotypic resistance (in blue). Phase III involves fitting regression models (discussed in the next section) to isolate metagenomic features that better correlate with variations of temperature and humidity (in red). AUC, area under the curve.

Fig. 3. ML performance and feature selection from correlations between gut microbial species, resistome and antibiotic resistance in E. coli.

a, Performance of the ML-powered predictive functions of E. coli resistance to specific antibiotics (ML technology: extra tree classifier; see Methods). Performance indicators (AUC, accuracy and precision) were computed as the average of 30 iterations of nested cross-validation (see Methods). See Supplementary Fig. 2 for performance indicator sensitivity, specificity and Cohen’s kappa score. The violin plots show the distribution of the data, with each data point representing one antibiotic model. Inside each violin plot is a box plot, with the box showing the interquartile range (IQR), the whiskers showing the rest of the distribution as a proportion of 1.5 x IQR and the white circle representing the median value. b, Counts of metagenomic features (ARGs and microbial species) found as the strongest predictors of E. coli resistance/susceptibility profiles to each antibiotic. c, Undirected graph showing the strongest predictors (metagenomic features in the chicken gut) for each antibiotic model. The edges of the graph link ARG or bacteria species nodes (predictor variables) to the antibiotic model in which they were found to be predictive. Both the ARG and antibiotic model nodes are colour coded according to the antibiotic class that the antibiotic/ARG is known to be associated with. The ML models were run for the following antibiotics: amoxicillin–clavulanic acid (AMC), amikacin (AMI), aztreonam (AZM), ceftazidime (CAZ), ceftazidime–clavulanic acid (CAZ-C), cefotaxime (CTX), cefotaxime–clavulanic acid (CTX-C), cefoxitin (CFX), chloramphenicol (CHL), cefepime (FEP), gentamycin (GEN), kanamycin (KAN), minocycline (MIN), nalidixic acid (NAL), streptomycin (STR), sulfafurazole (SUL) and trimethoprim–sulfamethoxazole (SXT). MDR, multidrug resistant.

Fig. 4. Gut features identified as predictors of E. coli resistance.

a,b, Microbial species and ARGs correlated with humidity (a) and temperature (b). Microbial species and ARGs are correlated with humidity or temperature, and also with each other, indicating that the ARGs are likely to be present in the species. Features were considered correlated if the slopes of the linear regression lines were significantly different from zero (P < 0.05 using a two-sided t-test). Nodes indicate ARGs or microbial species; edges connect species to ARGs likely present in the species. ARG nodes are colour-coded according to the antibiotic class known to be associated with the ARG; microbial species nodes are shown in grey.

Fig. 5. Gene structure and evolutionary analysis of the potentially mobile ARG pattern ISAba125–NDM-1.

Bayesian evolutionary phylogenetic tree reconstructing the phylogeny of contigs containing the clinically important ARG bla_NDM-1 and MGE ISAba125. Sample IDs (for example, LNPCJFT2-17) are given under the phylogenetic tree. The source type and location of the samples are indicated by coloured strips. The gene structure of each sample is shown below the tree with MGEs coloured blue, ARGs coloured green and other genes coloured yellow. The ARG ble is co-located with bla_NDM-1 in all contigs.