A reference gene catalogue of chicken gut antibiotic resistomes

Background: The chicken gut microbiota, as a reservoir of antibiotic resistance genes (ARGs), poses a high risk to humans and animals worldwide. Yet a comprehensive exploration of the chicken gut antibiotic resistomes remains incomplete. Results: In this study, we established the largest chicken gut resistance gene catalogue to date through metagenomic analysis of 629 chicken gut samples. We found significantly higher abundance of ARGs in the Chinese chicken gut than that in the Europe. tetX , mcr , and bla NDM, the genes resistant to antibiotics of last resort for human and animal health, were frequently detected in the Chinese chicken gut. The abundance of ARGs was linearly correlated with that of mobile genetic elements (MGEs). The host-tracking analysis identified Escherichia , Enterococcus , Staphylococcus , Klebsiella , and Lactobacillus as the major ARG hosts. Especially, Lactobacillus , an intestinal probiotic, carried multiple drug resistance genes, and was proportional to IS Lhe63 , highlighting its potential risk in agricultural production processes. Conclusions: We first established a reference gene catalogue of chicken gut antibiotic resistomes. Our study help to improve the knowledge and understanding of chicken antibiotic resistomes for knowledge-based sustainable chicken meat production.


Introduction
Antibiotic resistance has been recognized as a significant global threat to public health and food safety by the World Health Organization [1]. Currently, hundreds of thousands of people die each year from infections caused by antibiotic-resistant bacteria. Irrational use of antibiotics is thought to be the main driver for the development of antibiotic resistance [2]. Antibiotic resistance genes (ARGs) are widely distributed in surface water, sewage treatment plant effluent, soil, and animal feces [3]. In livestock farming environments, ARGs and drug-resistant bacteria can be transmitted to humans through the food chain, water, or air. Besides, the horizontal gene transfer (HGT) process makes it possible for bacteria from various ecological niches to share genes, thereby playing a vital role in the dissemination and evolution of ARGs. Undoubtedly, the gut atmosphere offers a strong communication site. A significant number of ARGs have been identified in the poultry gut [4].
As the world's population grows, the demand for animal protein including poultry meat, as a high-quality protein, is increasing day by day. Broiler chickens are considered one of the world's main meat suppliers. Global poultry meat consumption increased from 11 kg per capita in 2000 to 14,4 kg per capita in 2011 and is projected to reach 17,2 kg per capita in 2030 [5]. The production of poultry meat will reach 134.5 million tons in 2023, according to a report submitted by the Organization for Economic Cooperation and Development and the Food and Agriculture Organization [6]. The modernization of livestock farming is inseparable from the use of antibiotics, and it is estimated that the global consumption of antimicrobials will increase by 67% between 2010 and 2030 [7]. The continued use of antibiotics would alter the microbial structure of the animal gut and cause an increase in abundance of antibiotic resistance genes in the animal gut. This anthropogenic activity promotes the emergence and spread of antibiotic resistance genes including tetX, mcr, and blaNDM, the genes resistant to antibiotics of last resort for human and animal health.
Tigecycline is one of the last resorts for the treatment of complex infections caused by multi-resistant Gram-negative and Gram-positive bacteria. In 2019, two variants of plasmid-encoded tetX, named tetX3 and tetX4, were identified in Enterobacter and Acinetobacter isolated from animal and human sources. Plasmids carrying tetX3 and tetX4 are already widely distributed in human pathogens and are transferable and stable [8]. Polymyxin is an important lipopeptide antibiotic that serves as the last line of defense against multidrug-resistant gram-negative bacterial infections [9]. In November 2015, the first plasmid-mediated polymyxin resistance gene, mcr-1, was discovered in China [10]. mcr-1 and its variants have been identified in more than 50 countries to date. The mcr-1 gene is carried by multiple plasmids, including IncX4, IncHI2, and IncI2, which strengthens the risk of further mcr-1 transmission worldwide. The blaNDM gene encodes a carbapenemase that hydrolyzes almost all β-lactam antibiotics. Carbapenem-resistant Enterobacteriaceae that produce the blaNDM gene are found worldwide and have become a major public health threat.
Due to the diversity and abundance of microorganisms in the gut microecology, the human and animal gut is an important site for the increase and spread of ARGs, with a large number of ARGs settling and spreading to the gut flora through host cells.
Currently, studies of gut resistance groups have been conducted in a variety of animals, such as humans, pigs, and chickens [11], which have deepened the understanding of the distribution, evolution, and transmission of ARGs in the gut flora.
However, in general, there is still no standard process for metagenomic analysis.
In this study, we explored the comprehensive profile of ARGs in chicken gut by using large-scale metagenomic datasets. We established the largest chicken gut resistance gene catalogue to date through metagenomic analysis of 629 chicken gut samples. We analyzed the diversity, abundance of ARGs and ARG hosts from multiple perspectives. Besides, we reconstructed 4062 high-quality metagenome-assembled genomes (MAGs) using a binning approach, most of which belong to species lacking culture representation. Combining these individual efforts and establishing a unified non-redundant dataset of the chicken gut genome is essential to advance future microbiome studies. We think these data will lay the foundation for future studies on the structure of chicken gut microbiota.

Metagenomic sample collection
The 36 gut samples of chicken were collected from our previous study [12].
There were 60 chicken gut samples from the current laboratory study. A total of 533 chicken gut samples from 11 countries were randomly selected from NCBI, including Denmark (n = 20) and Bulgaria (n = 19) and USA (n = 11). Furthermore, all the above metagenomic datasets that were sequenced using the Illumina HiSeq platform with PE150 (paired-end sequencing 150 × 2) were downloaded to avoid the fluctuations caused by different sequencing strategies. Additional sample descriptions can be found in Supporting Information Fig. 1 and Table S1.

Quality control and metagenome assembly
First, raw reads with average quality scores < 20 (Q20) or length < 20 bp (L20) were removed using Sickle [13] (https://github.com/najoshi/sickle). After that, the metagenomic sequencing data of each sample were de novo assembled with default k-mer size using the CLC Genomics Workbench (Version 10.0.1, CLC Bio, Aarhus, Denmark), and the assembled contigs longer than 500 bp were reserved. We finally obtained 35672711 contigs from 629 metagenomic samples. The detailed information is summarized in Table S2.

Identification of ARG-like ORFs and the mobile genetic elements (MGEs)
The open reading frames (ORFs) of the contigs were predicted using Prodigal v2.6.3 with a metamodel [14]. The ARG-like ORFs were identified against the resfinder database [15] using BLASTN with an E-value of ≤ 1×10 -5 , and the results that identity ≥ 80% and over 70% query coverage were identified as ARG-like ORFs [16]. The coverage of ARG-like ORFs was calculated by mapping metagenomic reads to the contigs with a minimum length coverage of 95% at 95% similarity and contigs minimum length ≥ 500bp using the CLC Genomics Workbench [16]. The number of cells in each metagenome sample was calculated by ARGs-OAP v2.0 [17] ( Table S1 for complete sample cell count information). To compare the coverage of each different sample, the ARG-like ORFs coverage was normalized by the number of cells in each sample (copies/cell). The calculation formula was as follows: where N is the number of reads mapped to ARG-like ORFs, L is the sequence length of the target ARG-like ORFs, n is the number of ARG-like ORFs, and 150 is the length of Illumina sequencing reads [12,16], and the C means the cell number of per sample. The final abundance calculation method was the same as the above equation remains the same.

Phylogenetic analyses
All phylogenetic trees were inferred using MEGA X (Version 10.

Identification of bacterial host of ARG-carrying contigs
The predicted protein sequences of ORFs from the contigs that carried ARG-like ORFs were annotated using BLASTP against non-redundant protein NCBI databases

Statistical analysis
The heatmap, correlation heatmap, and variance analysis were generated using the pheatmap, ggcorrplo, and limma packages in the R package (v4.0.2, https://www.r-project.org/), respectively. The alpha diversity of the gut microbiota, including the number of species and the Shannon index, was calculated by spaa and vagen in the R package (v4.0.2). MetaPhlan2 (v2.7.5)[31] was used to identify species and to determine their relative abundances across samples. The network was explored and visualized using the interactive platform of Gephi (v0.9.2). The drawing of simple graphs was done by GraphPad Prism 8. Statistical analyses were performed using SPSS software (PASW Statistics 22.0). The Pearson correlation coefficient was used to show the correlation between the abundance of ARGs and the abundance of MGEs.
A two-tailed Wilcoxon test (two-by-two comparison) was used to compare Chinese and European chickens. A P-value < 0.05 was considered a statistically significant difference. The Procrustes transformation was performed using two Bray-Curtis distance plots (PCoA) as input based on microbial community matrices and ARGs at the subtype level. R (v4.0.2), Python (v3.7.3) programs were used to handle the text processing, information extraction, and data statistics for this study.
Depending on the geographical location of the samples, we could divide them into two groups, China and Europe (Due to the small number of samples in the USA, there were no group discussions). Overall, the average abundance of ARG types in the chicken gut was 3.62 copies/cell. The tetracycline (1.24 copies/cell), aminoglycoside (0.70 copies/cell), and macrolide (0.69 copies/cell) resistance genes were the predominant ARG types. (Fig. 2b). The abundance of ARG types in China was significantly different from that in Europe (P < 0.0001). The average abundance of ARGs was about three times higher in China than that in Europe (Fig. 2c). Human activities (different strategies for antibiotic use in different geographic areas) may be the main driver of the different distribution of ARGs.
The similarity of ARG subtypes compositions was calculated using principal component analysis (PCA). Chicken gut samples were well clustered in Europe, whereas in China they were more dispersed and well surrounded by European samples, suggesting that the Chinese region has a more diverse ARG composition than that in Europe (Fig. 2a). Each sample carried a different number of species of ARGs, so we further compared the diversity of ARGs between China and Europe.
Interestingly, the results showed the presence of more diverse resistance gene types in Europe ( Fig. 2d and Fig. 2e). This may be related to the fact that the samples collected in Europe were spread over nine countries.
Of the 387 ARG subtypes detected in the chicken gut, the top five abundant ARG subtypes were considered representative ARGs, tetW, ermB, tetA, tetL, and sul2.
Interestingly, we found that more than half of these five dominant ARG subtypes were tetracycline resistance genes, of which tetW accounted for 18.25% of the overall.
Concerning tetX genes, tetX10 and tetX2 were more prevalent in the chicken gut and most were present on plasmids. Most of tetX belonged to the Bacteroides and Butyricimonas. IS4351, IS1187, IS612, IS614, and ISBfl1 may have an important role in the transfer of tetX. Most of mcr-1 were present on plasmids that were attributed to Escherichia coli, Escherichia albertii, and Klebsiella pneumonia. As well, ISApl1 was the predominant transposon for mcr-1 transfer. Most of blaNDM were also present on plasmids carried by a variety of pathogenic bacteria such as E. coli, Salmonella enterica, and Proteus mirabilis. ISABa125, IS5D, and IS26 were potential for their transfer vectors (Fig. 3).

Bacterial community structure
Overall, the Firmicutes (59.13%), Proteobacteria (25.38%), Bacteroidetes Escherichia were the two genera with the most significant differences between the Chinese and European groups, respectively (Fig. 4c). The correlation analysis revealed that the abundance of Lactobacillus and Escherichia were negatively correlated (ρ = -0.503,P < 0.01, Spearman's correlation coefficient, two-tailed test).
Previous studies have shown that bacterial community structure determines the distribution of resistance genes without considering HGT, and a significant correlation between bacterial composition and ARG profiles could observe [32,33]. Visualized by the Procrustes analyses, ARGs and bacterial composition in the chicken gut showed highly significant goodness of fit by ANOVA (M 2 = 0.7888, P < 0.001) (Fig.   S1).

The host of ARGs in chicken gut samples
After filtering the contigs with lengths less than 500bp, we obtained a total of 35672711 contigs from 629 samples, of which 24318 contigs were identified as ARG-carrying contigs (ACCs) ( ACCs. The phylum with the highest percentage of these 6409 contigs with coexisting ARGs and MGEs was Proteobacteria (62.16%), while at the level of the genus, the highest percentage was Escherichia (12.65%), followed by Klebsiella(6.54%). Also in terms of the types of drug resistance included, the highest percentage was aminoglycosides (26.56%), followed by phenicol, β-lactam, and sulfonamide with 16.54%, 13.90%, and 13.48%, respectively. We further explored the correlation between intI1 and ARGs. intI1 was found to be correlated with aph(6)-Id, blaOXA-10, blaTEM-1B, mphA, floR, qnrS1, sul1, sul2, sul3, and tetA, which are drug resistance gene subtypes significant positive correlations (Table S4).
We found a linear correlation between the abundance of MGEs and the abundance of ARGs (y = 3.15+1.35 x, R 2 = 0.3) (Fig. 6a), and also the correlation heat map produced by ggcorrplot package showed a correlation between ARGs and MGEs (r = 0.56, P < 0.05) (Fig. 6c), indicating that the presence of MGEs plays a transfer of ARGs an important role. We found the highest abundance of ISLhe63 among the IS and ISLHe63 showed a significant correlation with Lactobacillus(ρ = 0.621, P < 0.05), which is consistent with previous studies[36].
Using ARGs structure as the response variable, the top 10 genus in terms of abundance and 3 types of MGEs(Integron, IS, Plasmid) as explanatory variables for RDA, after forward model selection, provided 12 explanatory variables, all of which were significantly positively correlated except for Bacteroides (Fig. S2a). VPA was used to determine the co-occurrence of MGEs and genus on ARGs structure. The results showed that the most important influences on the structure of ARGs were MGEs, with individual explanatory rates of 12.1% and 6.6%, respectively (Fig. S2b).

MAGs from the chicken gut
On each of the 629 samples, we carried out a single-sample metagenomic assembly. A total of 25598 MAGs were created from binning. With a total of 4062 high-quality genomes with completeness ≥ 80% and contamination ≤ 10%, we were left. Of these 1823 had completeness > 90% and contamination < 5%, 938 were > 95% complete with < 5% contamination, and 163 MAGs were > 97% complete with 0% contamination. The distribution of these MAGs between the 629 samples can be found in Table S5.  Due to the limitations of the taxonomically annotated reference sequences, some MAGs failed to classify to the genus level. This limitation is common in the taxonomic identification of genomes generated from metagenomic data [29]. Despite this, this study still revealed that 156 bacteria genera were the hosts of 140 ARG subtypes. Therefore, this study greatly expands the knowledge of the diversity of ARBs in the chicken gut.

Discussion
Our study analyzed antibiotic resistance based on metagenomic tools within 629 chicken gut samples, demonstrating the correlation between ARGs and MGEs and revealing ARG hosts from multiple perspectives. The microbiome in the chicken gut microbiota constitutes a large reservoir of ARGs. Tetracyclines, aminoglycosides, and macrolides were the most commonly observed ARG forms in this study, whereas certain resistance genes are often present in particular habitats, such as urban We found that the total level of acquired antimicrobials correlates with the overall livestock antimicrobial use. The average ARG abundance in the gut of chickens in China was three times higher than that in Europe. According to a survey report by the World Organization for Animal Health (OIE) in February 2019, Asia, the Far East, and Oceania had the highest average use of antimicrobials in animals in the world, amounting to 296.75 mg/kg, which is 3.74 times higher than in Europe. This value is close to the average abundance ratio of ARGs between China and Europe above. The total ban of all antibiotic growth promoters in the EU member states was in 2006 while it in China was in 2020. Different strategies of human antibiotic use in agricultural practices may be responsible for this phenomenon.
ARGs were frequently detected in the chicken gut, such as tetX, mcr, and blaNDM, which may be due to the selection posed by prolonged exposure of gut microbes to antibiotics. Although tigecycline was not applied to livestock animal farming, tetX was still found in the gut environment of chickens in this study, presumably due to the selection pressure of tetracyclines. The latest appearance from animals and humans of plasmid-mediated high levels of tigecycline resistance genes tetX3 and tetX4 has triggered the much broader public alarm [8], rendering it important for robust surveillance of tetX genes to be improved. In November 2015, the plasmid-mediated polymyxin resistance gene mcr-1 was first identified, subsequently, multiple Enterobacteriaceae members with different host origins were found carrying the resistance gene in more than 50 countries per continent, suggesting that mcr-1 was spread globally [41]. And it has been demonstrated that the transposon part Tn6330 We constructed 4062 high-quality MAGs from 629 chicken gut samples, greatly extending the enrichment of previous chicken gut MAGs [29,52]  However, it has also been shown that multiple ARGs are isolated from Lactobacillus strains and some of them may even transfer their intrinsic ARGs to other lactic acid bacteria or pathogens through horizontal gene transfer [64][65][66][67]. In our study, Lactobacillus was also found to carry a variety of ARGs, especially macrolides and tetracyclines. This suggests that Lactobacillus can potentially act as a reservoir of ARGs and may play an active role in the transfer of ARGs to humans through the food chain. We identified 4 archaea in the 4062 high-quality MAGs. It has been shown that Methanobrevibacter and Methanocorpusculum were positively correlated with all of the detected ARGs [68]. Although we did not identify any ARGs in the 4 archaeal strains, it still provides a new idea for a potential link between archaea and ARGs.

Conclusion
We established the largest chicken gut resistance gene catalogue to date through metagenomic analysis of 629 chicken gut samples and revealed the diversity and abundance of ARGs in the chicken gut. The results explained that the differences in the abundance levels of ARGs (e.g., tetX, mcr, and blaNDM) between different regions in China and Europe were completely positively correlated with anthropogenic factors.
The widespread presence of MGEs further drives the spread of ARGs through horizontal transmission. Lactobacillus, as gut probiotics, also carried multiple drug resistance genes (especially macrolides and tetracyclines), emphasizing their potential risk. Besides, the number of chicken gut microbial genomes in public databases was enriched by reconstructing a large number of MAGs, but a whole-genome sequencing approach after isolation and culture is still needed to validate the description. Our study helps to improve the knowledge and understanding of chicken gut flora and promote knowledge-based sustainable chicken production.

Consent for publication
Not applicable.

Data availability
The dataset supporting the conclusions of this article is included within its supplementary file.

Competing interests
The author declares that have no competing interests.

Funding
This work was supported by the National Natural Science Foundation of China (no.   Table S1. Sample information.         represents the log2-fold change and the x-axis represents the genus name.

Fig. 5
The taxonomy of ACCs (in genus level) and the percentages of ARG types these contigs carried. For example, Escherichia (6.3%) represents that 6.3% of ACCs were annotated as Escherichia. b Bar chart shows the percentages of ARG types that were carried by the annotated ACCs.For example, 44.1% of the ARG-carrying contigs originating from Escherichia carried aminoglycoside resistance genes.   Figure 1 Distribution map of chicken gutl samples and diversity of ARGs, pie chart showing the pro le of ARG abundance (top 3 ARG types). Europe(n=178) includes Netherlands(n=20), Germany(n=19), France(n=20), Spain(n=20), Belgium(n=20), Italy(n=20), Poland(n=20), Denmark(n= 20) and Bulgaria (n=19). ARGs heatmap and phylogenetic tree of tetX, mcr, and blaNDM detected in the chicken gut.  The taxonomy of ACCs (in genus level) and the percentages of ARG types these contigs carried. For example, Escherichia (6.3%) represents that 6.3% of ACCs were annotated as Escherichia. b Bar chart shows the percentages of ARG types that were carried by the annotated ACCs.For example, 44.1% of the ARG-carrying contigs originating from Escherichia carried aminoglycoside resistance genes.  The network analysis revealing the co-occurrence patterns between major MAGs(in genus level) and ARG subtypes. The nodes were colored according to the modularity class. Modularity index 0.68. The size of each node was proportional to the number of connections, that is, the average weighted.