Whole Resistome Analysis in Campylobacter jejuni and C. coli Genomes Available in Public Repositories

Campylobacter spp. are the most frequent agent of human gastroenteritis worldwide, and the spread of multidrug-resistant strains makes the clinical treatment difficult. The current study presents the resistome analysis of 39,798 Campylobacter jejuni and 11,920 Campylobacter coli genomes available in public repositories. Determinants of resistance to β-lactams (Be) and tetracyclines (Te) were the most frequent for both species, with resistance to quinolones (Qu) as the third most important on C. jejuni and to aminoglycosides (Am) on C. coli. Moreover, resistance to Te, Qu, and Am was frequently found in co-occurrence with resistance to other antibiotic families. Geographical differences on clonal complexes distribution were found for C. jejuni and on resistome genotypes for both C. jejuni and C. coli species. Attending to the resistome patterns by isolation source, three main clusters of genomes were found on C. jejuni genomes at antimicrobial resistance gene level. The first cluster was formed by genomes from human, food production animals (e.g., sheep, cow, and chicken), and food (e.g., dairy products) isolates. The higher incidence of tet(O), associated with tetracycline resistance, and the gyrA (T86I) single-nucleotide polymorphism (SNP), associated with quinolone resistance, among genomes from this cluster could be due to the intense use of these antibiotics in veterinary and human clinical settings. Similarly, a high incidence of tet(O) genes of C. coli genomes from pig, cow, and turkey was found. Moreover, the cluster based on resistome patterns formed by C. jejuni and C. coli genomes of human, turkey, and chicken origin is in agreement with previous observations reporting chicken or poultry-related environments as the main source of human campylobacteriosis infections. Most clonal complexes (CCs) associated with chicken host specialization (e.g., ST-354, ST-573, ST-464, and ST-446) were the CCs with the highest prevalence of determinants of resistance to Be, Qu, and Te. Finally, a clear trend toward an increase in the occurrence of Te and Qu resistance determinants on C. jejuni, linked to the spread of the co-occurrence of the blaOXA–61 and tet(O)-tet(O/W/O) genes and the gyrA (T86I) SNP, was found from 2001 to date in Europe.


INTRODUCTION Thermophilic
Campylobacter species, and specially Campylobacter jejuni and Campylobacter coli, are the most frequent agents of human gastroenteritis both in industrialized countries (Kaakoush et al., 2015;European Food Safety Authority, 2018) and in low-and middle-income countries (Fischer Walker et al., 2013;Liu et al., 2015), with C. jejuni accounting for up to 80% of all human Campylobacter gastroenteritis cases (Bronnec et al., 2016), with estimations of over 800,000 cases per year between 2000 and 2008 (Scallan et al., 2011). Moreover, C. jejuni is the most frequent foodborne bacterial pathogen globally, with developing countries being the most affected (Yahara et al., 2017).
Campylobacter jejuni and C. coli are present in the digestive tract of different animals, including birds, cattle, sheep, pigs, and pets (Thépault et al., 2017), although chicken is recognized as the main source of human infection in most countries (Domingues et al., 2012;Newell et al., 2017;Cody et al., 2019). Despite consumption of contaminated raw or under-cooked poultry meat being the main source of human infection (Domingues et al., 2012), other sources of infection include contaminated water, unpasteurized/incompletely pasteurized milk, or the direct contact with animals and the environment (Lévesque et al., 2008;de Perio et al., 2013;Levallois et al., 2014;Fernandes et al., 2015;Mossong et al., 2016).
Antimicrobial-based therapy for campylobacteriosis is indicated only in patients with severe and invasive gastroenteritis. Macrolides (mainly erythromycin), fluoroquinolones (mainly ciprofloxacin), and tetracyclines are the most commonly used and effective antibiotics (Aarestrup et al., 2008;Silva et al., 2011). Unfortunately, these antimicrobial agents have analogs widely employed in veterinary settings, and their overuse or misuse have favored the emergence and spread of resistant Campylobacter isolates (Van Boeckel et al., 2015;Asuming-Bediako et al., 2019). Furthermore, the continuous increase in resistance to certain antimicrobial classes has been accompanied by a rise in the occurrence of multidrug-resistant (MDR) C. jejuni and C. coli from farm to fork (Mourkas et al., 2019).
Whole genome sequencing (WGS) is a powerful tool that offers high-resolution sub-typing of Campylobacter isolates, which is of value for outbreak investigations (Clark et al., 2016;Lahti et al., 2017). Apart from multilocus sequence typing (MLST) schemes, based on the analysis of the sequence of several house-keeping genes, WGS also allows to predict relevant phenotypes, such as antibiotic resistance, in Campylobacter (Zhao et al., 2016). Moreover, recent studies have found high concordance between the presence of antimicrobial resistance determinants in WGS data and phenotypic resistance among Campylobacter isolates from Latvia (Meistere et al., 2019), England and Wales (Painset et al., 2020), and Poland (Wysok et al., 2020), supporting the use of WGS data as a good predictor for phenotypical antimicrobial resistance.
Taking into account the increasing relevance of C. jejuni and C. coli as foodborne pathogens and the availability of a high amount of their genomes in different public repositories, with the corresponding information on the year, country, and source of isolation of the sequenced strains, the objectives of the current study are as follows: (i) to study the host association of the most abundant clonal complexes (CCs) of C. jejuni and the most abundant sequence types (ST) within the main CCs on C. jejuni and C. coli; (ii) to compare the prevalence of antimicrobial resistance (AR) factors within C. jejuni and C. coli genomes; (iii) to perform an exhaustive resistome genotypic analysis and link AR patterns to country/continent, source of isolation, and C. jejuni CCs; and (iv) to assess temporal changes in such AR patterns along the last 20 years.

Downloading Genomes and Metadata for Each Genome
Metadata associated to each genome was obtained by downloading the table obtained after using "Campylobacter jejuni" and "Campylobacter coli" as keywords (i) on PATRIC Taxon View website 1 ; (ii) through the Export dataset section on the PubMLST-Campylobacter jejuni/coli website 2 , where all isolates were selected; and (iii) by using 02.NCBI_Download_metadata.R script 3 for the NCBI repository. This script basically downloads the information obtained from https://www.ncbi.nlm.nih.gov/biosample/ BIOSAMPLE?report=full&format=text, where BIOSAMPLE is substituted by the NCBI Biosample code for each genome analyzed, with these genomes being those of the C. jejuni and C. coli strains compiled in the prokaryotes.txt file in ftp://ftp.ncbi.nlm.nih.gov/genomes/GENOME_REPORTS/. The three metadata files obtained (for the NCBI, PATRIC, and PubMLST repositories) were manually inspected and adapted to merge them in a unique file. The collection_date, host_name, country, and continent columns were manually inspected to homogenize the collection date (to year of collection), isolation source, and location format (to country and continent).
C. jejuni and C. coli genomes publicly available at NCBI, 4 PATRIC, 5 and PubMLST 6 repositories were downloaded on February 2021 using the download_genomes pipeline made for such task (see footnote 3). Briefly, the lists of available C. jejuni and C. coli genomes were constructed from metadata files for the PATRIC and PubMLST repositories and from the prokaryotes.txt file 7 for the NCBI repository. Those genomes that were present on more than one repository were discarded from the NCBI or PubMLST dataset using 03.ena2ncbi.rb and 04.remove_repetition scripts (see footnote 3), and then, all selected genomes were downloaded using the corresponding script for each database from download_genomes github repository (see footnote 3).

Genome Analyses
All genomes were analyzed through MLST by using the mlst github pipeline (Seemann), which employs data from the PubMLST.org website (Jolley et al., 2018). Antimicrobial gene detection and detection of single-nucleotide polymorphisms (SNPs) conferring antimicrobial resistance were performed using the ResFinder database version updated at March 12, 2021 (Zankari et al., 2012) and PointFinder database updated at February 1, 2021 (Zankari et al., 2017), respectively. Both analyses, together with MLST, were performed for each genome by using the staramr pipeline (phac-nml/staramr) and 08.staramr_auto.rb script (see footnote 3) to automatize the analyses for the dozens of thousands genomes analyzed.

Statistical Analysis
The R-packages tidyr and dplyr were used to organize the data matrix, while the R-packages ggplot2 and pheatmap were employed for the representation of line plots, boxplots, and heatmap plots using R version 3.6.1 (R Core Team, 2019) and t-test to determine statistical differences. Heatmap plots were done using the percentage of genomes carrying each AR determinant (gene or SNP) that belonged to the same isolation source, clonal complex, or isolation year. The Euclidean distance and the complete (= UPGMA) clustering method were employed for the clustering of columns and rows.
Only those isolation sources with more than 50 genomes, excluding the unknown group, were selected for obtaining heatmap plots of antibiotic resistance genes (ARGs) or antibiotic families versus isolation source to avoid the introduction of biases due to the low number of genomes. Similarly, European genomes were uniquely selected for temporal variation analyses in the resistome, due to the fact that the highest number of genomes belonged to Europe and North America, but an important proportion of North-American genomes had missing information in relation to the year of isolation. Moreover, only European genomes from isolates obtained in those years when at least 100 genomes were available (i.e., years between 1997 and 2018, with the exception of 1999, 2000, and 2002) were selected for this task.

Global Overview
A total of 39,798 C. jejuni and 11,920 C. coli genomes were analyzed. The C. jejuni genomes belonged to 41 CCs, with ST-21 CC being the most abundant (23.8% of the analyzed genomes) followed by  and 7.8% of the analyzed genomes, respectively), while C. coli genomes belonged to 21 CCs, with ST-828 CC as the main CC (84.8% of analyzed genomes). Remarkably, 13.7% of the C. jejuni and 13.43% of C. coli genomes were not assigned to any known CC. The main isolation source was humans (37.1% of C. jejuni genomes and 14.8% of C. coli genomes), followed by chicken, cows, and sheep for C. jejuni (10, 3, and 1.2% of genomes, respectively) and chicken, pigs, cows, and sheep for C. coli (9, 5.5, 2.3, and 1.3% of genomes, respectively). Unfortunately, 45.6 and 60.6% of C. jejuni and C. coli genomes, respectively, had missing information on the isolation source. A total of 18,812 C. jejuni (47.3%) and 7,959 C. coli (66.8%) genomes were from North-American isolates, while 18,787 C. jejuni (47.2%) and 3,328 C. coli (27.9%) genomes belonged to European isolates. The genome codes, collection year, continent, country, host and MLST information for all the analyzed C. jejuni and C. coli genomes were presented to Supplementary Tables 1, 2, respectively.
The most abundant AR factors were associated with resistance to β-lactams and tetracyclines for both species, although β-lactam resistance was most prevalent on C. jejuni (81.6 vs. 65.0%). AR factors against aminoglycosides and antibiotics of the macrolides-lincosamides-streptogramins-pleuromutilins group (MLSP) were most abundant on C. coli genomes (34.8 and 2.9%, respectively) than on C. jejuni (6.7 and 0.5%, respectively), while resistance to quinolones was much more frequent on C. jejuni (28.7%) then on C. coli (0.6%) genomes (Figure 1). Moreover, up to 7.8 and 23.3% of C. jejuni and C. coli genomes, respectively, did not harbor any known AR factor (Figure 1). Remarkably, 39.9 and 21% of C. jejuni and C. coli genomes, respectively, had only AR genes or SNPs associated with resistance to only β-lactams, while 41.7 and 44.0% of the genomes harbored determinants of resistance to β-lactams combined with determinants of resistance to other antibiotic families. Similar observations were made for resistance to tetracyclines, which were seen as unique resistance factors on 4.6 and 5.1% of C. jejuni and C. coli genomes, respectively, while 38.6 and 41.6% of genomes harbored resistance to tetracyclines combined with resistance to other antibiotic families (Figure 1). It is worth highlighting the high level of co-occurrence of resistance to aminoglycosides, β-lactams, and tetracyclines on C. coli genomes (20.2%) and to β-lactams, quinolones, and tetracyclines on C. jejuni genomes (15.4%) (Figure 1).
A total of 120 and 68 different AR genes and SNPs were detected for C. jejuni and C. coli genomes, respectively. bla OXA−61 and tet(O) genes were the most extended AR determinants on both C. jejuni and C. coli genomes (they were carried by 69.5 and 36.4% of C. jejuni and 35.8 and 43.6% of C. coli genomes, respectively). The gyrA (T86I) SNP was the third most prevalent AR factor on C. jejuni genomes (28.5%), while aph(3 )-III, bla OXA−489 , aadE-Cc, and bla OXA−193 were present in more than 10% of the C. coli genomes (Figure 1).
Moreover, while 33.48% of C. jejuni genomes carried only the bla OXA−61 gene, 10.7% simultaneously harbored bla OXA−61 and tet(O), and 9.8% had gyrA (T86I), bla OXA−61 , and tet(O) in combination, with these three ARG genotype patterns being the most abundant on C. jejuni genomes ( Table 1). The presence of only the bla OXA−61 gene was observed on 10.4% of the C. coli genomes, with it being the main genotype pattern found, followed by only the bla OXA−489 gene and bla OXA−61 together with tet(O), with 6.9% and 6.6% of prevalence, respectively ( Table 1).

Geographical Distribution of CCs and ARGs-SNPs
The first global analysis of C. jejuni CC distribution showed differences by continent. The ST-21 CC was the most abundant FIGURE 1 | Global antibiotic resistance gene (ARG)-single-nucleotide polymorphism (SNP) distribution on C. jejuni and C. coli genomes. Antimicrobial resistance families, family co-occurrence, and most abundant ARGs or SNPs among the C. jejuni and C. coli genomes analyzed. Am, aminoglycosides; Be, beta-lactams; Qu, quinolones; Te, tetracyclines. Gene names are colored according to the corresponding antimicrobial family [blue: beta-lactams, orange: tetracyclines, red: aminoglycosides, purple: quinolones, and green: macrolides-lincosamides-streptogramins-pleuromutilins group (MLSP)]. Hyphen (-) indicates no ARGs or SNPs detected.
in Asia, Europe, and North America (26.9, 25.0, and 23.0% of genomes, respectively), while ST-353 CC was predominant in South America and Africa (23.7 and 15.6%, respectively) and ST-354 CC was dominant in Oceania (60.3%) (Figure 1). These patterns are not observed in all the countries from the same continent. For example, ST-353 CC was abundant in the United States and Brazil (12.8 and 33.5% of genomes, respectively), but it was less represented on Canadian and Peruvian genomes (2.2 and 7.0%, respectively), with Canada having a high prevalence of ST-45 CC (18.8%) (Figure 2). The main differences among countries within the same continent were found in Europe, with 93.3% of ST-21 CC genomes in Denmark; 34.4% of ST-45CC genomes in Finland; and 26.7, 25, and 39.5% of ST-21 CC genomes in France, United Kingdom, and Israel, as the main CCs (Figure 2).
The bla OXA−61 gene had the highest prevalence among C. jejuni genomes in all continents (from 36.5 to 92.4% of total genomes), except for South America, which had the gyrA (T86I) SNP as the main AR factor (39.1%) ( Figure 3A). Oceania had the highest prevalence of tet(O), followed by North America and Asia (60, 46.9, and 42.7%, respectively), while the gyrA T86I SNP had the highest abundance on genomes obtained in Asia and Oceania (73.9 and 59.2%, respectively), which cluster together and separated from those from other continents ( Figure 3A and Supplementary Figure 1A). There were some discrepancies at country level, where two main clusters were obtained, the second one containing genomes from Taiwan, Israel, New Zealand, Luxemburg, and China, countries with no geographical or even socio-economical relation, which harbored the highest prevalences of tet(O) gene and gyrA (T86I) SNP ( Figure 3B and Supplementary Figure 1B).
The tet(O) and bla OXA−61 genes were the most abundant on C. coli regardless of the continent, with prevalences ranging from 16.9 to 36.1% and 13.7-28.1% of genomes, respectively. The aadE-Cc and aph(3 )-III genes showed an important presence on European genomes (22.6 and 13.6%, respectively), while ant(6')-Ia had an important prevalence on North and South American genomes (11.8 and 15.6%, respectively) and the gyrA (T86I) SNP on Asian genomes (12.2%) ( Figure 3C). At country level, using those countries with at least 20 genomes, the cluster formed by genomes from China, Egypt, Switzerland, Peru, Canada, United States, and Vietnam had a higher prevalence of aph(3 )-III, bla OXA−61 , and tet(O) genes than the cluster formed by genomes from the United Kingdom, Sweden, Australia, India, and Chile ( Figure 3D and Supplementary Figure 1).

Host Specialization and ARG-SNP Distribution by Source of Isolation
The analysis of the occurrence of C. jejuni CCs in each of the isolation sources, excluding human, showed three clear clusters of CCs. Chicken specialist and generalist CCs formed the Only ARG genotype patterns with relative abundance > 1% are presented.
biggest group, which was clearly subdivided into two sub-clusters according to host range. Interestingly, chicken was the most frequent source found on the generalist sub-cluster ( Figure 4A). Previous studies (French et al., 2005;Griekspoor et al., 2010;Revez et al., 2011;Sheppard et al., 2014;Dearlove et al., 2016;Mourkas et al., 2020) have classified up to eight C. jejuni CCs as host generalist and 13 CCs as host specialist, including three CCs specialized in cattle infection, seven in chicken infection, and other three in wild bird infection or environmental persistence (  Figure 4A and Table 2). On the other hand, ST-283 CC, previously described as a specialist chicken CC, was not located within the chicken specialist sub-cluster (Figure 4A), so a change to generalist host range can be proposed according to the data here presented. Another important observed cluster was the cluster of cattle specialist CCs, which also included the ST-22 CC, previously described as generalist ( Figure 4A). In this case, the percentage of genomes assigned to the cow source among cattle specialist CCs, which ranged from 44.1 to 66.7% (excluding genomes isolated from human), was not as high as that of genomes assigned to chicken source for chicken specialist CCs, which was in all cases above 90%. In addition, cattle specialist CCs had similar percentages of abundance in the cow source than generalist CCs in the chicken source ( Figure 4A and Supplementary Figure 2A). Therefore, it is important to remark the high level of specialization of those CCs considered as chicken specialist, compared to other specialist groups. Moreover, a ST host specialization analysis performed with the most abundant CCs showed that ST-45 CC harbors STs  Table 3. The analysis of the occurrence of AR determinants in genomes from each of the isolation sources revealed three main clusters of both C. jejuni and C. coli genomes. The main cluster of C. jejuni genomes related to food production sources, which showed a higher prevalence of bla OXA−61 (from 44.7 to 100%) and of the gyrA (T86I) SNP (from 1.5 to 41.1%) than genomes from the second cluster obtained [from 0 to 12.84% for bla OXA−61 ; from 1.3 to 4.6% for the gyrA (T86I) SNP] (Figure 5A and Supplementary Figure 3A). This second cluster harbored genomes from isolates coming from host animals not linked to industrialized food production, such as jackdaws, crows, and other wild birds ( Figure 5B). The lower prevalence of bla OXA−61 and the higher diversity of other genes associated with resistance to β-lactams among the genomes from this second cluster indicate a lower selection pressure to keep this prevalent bla OXA−61 gene in non-livestock animal populations. The last cluster of C. jejuni genomes was formed by genomes from duck, other animals, and environmental waters and presented an intermediate prevalence of bla OXA−61 compared to the other two clusters obtained (Figure 5A and Supplementary Figure 3A). Similar clusters were found for C. coli genomes, with pig, cow, and turkey forming the cluster with the highest prevalence of tet(O) and bla OXA−61 , with tet(O) being more prevalent than on C. jejuni genomes while bla OXA−61 was less prevalent (Figure 5B and Supplementary Figure 3B). A second cluster of C. coli genomes was formed by genomes from dairy products, humans, chicken, and other animals and was characterized by an intermediate prevalence of tet(O) and bla OXA−61 genes, compared to the other clusters ( Figure 5B and Supplementary Figure 3B). Finally, the third cluster, with the lowest prevalence of AR determinants, was formed by genomes from environmental samples, sheep, and duck. It is important to highlight the high prevalence of aadE-Cc on genomes from pigs and bla OXA−489 on genomes from sheep ( Figure 5B).

ARG-SNP Distribution Among C. jejuni CCs
Up to four main clusters were obtained in the analyses performed at antibiotic family level. Cluster 1 harbored four CCs with the highest prevalence of determinants of resistance to tetracyclines and quinolones, while cluster 2, formed by 20 CCs, had the second main prevalence for tetracycline and quinolone resistance determinants ( Figure 6A and Supplementary Figure 4A). Cluster 3, formed by seven CCs, had the lowest prevalence of quinolone and tetracycline resistance determinants ( Figure 6A and Supplementary Figure 4A), and cluster 4, formed by six CCs, had the lowest abundance of β-lactam resistance determinants ( Figure 6A and Supplementary Figure 4A). At gene level, three clusters were found: cluster 1, harboring four CCs with a high prevalence of gyrA (T86I) SNP, bla OXA−61 , and tet(O); cluster 2, composed of 20 CCs with a high prevalence of bla OXA−61 ; and cluster 3, with 14 CCs showing the lowest prevalence of gyrA (T86I) SNP, bla OXA−61 , and tet(O) (Figure 6B and Supplementary Figure 4B). No clear association between genotype patterns at antibiotic family level and host specialist/generalist status was observed. As an example, cluster 1, characterized by a high prevalence of ARGs linked to resistance to tetracyclines and quinolones, was mainly formed by chicken specialist CCs, but other chicken specialist CCs were located across the other three remaining clusters (Figure 6A). A similar absence of a Genomes belonging to the human isolation source were removed before the calculation of the relative abundances, and CCs/STs with 20 or less genomes were removed to avoid biases due to the low number of genomes. The Host range classification of C. jejuni CCs is according to previous publications (see Table 1), except for those with highlighted edges, where edge color indicates previous classification while fill color indicates the new assignment, according to the results obtained in the current study.
clear association with the host range of each CC was observed at gene level.

Temporal Changes of AR Determinants on Europe
Focusing on temporal changes in the repertoire of AR determinants among the C. jejuni European genomes, there was a clear stability in the occurrence of β-lactam ARGs along time (with mean prevalence values per each 5-year period ranging from of 88.8 to 93.3%), while an important increase was observed in the prevalence of determinants of resistance to tetracyclines, from a mean value of 29.  Figure 5B), which were the main determinants of resistance to tetracyclines and quinolones, respectively. Furthermore, the increases along time in the tet(O/W/O) gene and the gyrA (T86I) SNP were observed in both human and chicken genomes (Figure 7C and Supplementary Figures 5C,D). Finally, the enhanced occurrence over time of resistance to tetracyclines and quinolones was mainly correlated with a significant increase in bla OXA−61 -tet(O)-gyrA (T86I) and bla OXA−61 -tet(O/W/O)-gyrA (T86I) multi-drug resistance genotype patterns on C. jejuni genomes obtained from human and chicken, which showed mean prevalence values of 10.5 and 6.1% in the 2000-2004 5-year period and 25.6 and 29.9% in the 2015-2018 period, respectively, representing a significant increase in co-occurrence of such genes ( Figure 7D and Supplementary Figures 5E,F).
Similar temporal analyses were performed on C. coli genomes, but no clear trends of AR determinant increase along time were detected. The observed slight increases in the prevalence of determinants of aminoglycoside and tetracycline resistance (Figure 8A), highly correlated with an increase of addE-Cc and tet(O) prevalence (Figure 8B), were associated with a differential prevalence on genomes from cows and pigs, with more C. coli genomes being sequenced on the last years (Figures 8C,D).

DISCUSSION
At a first global view, the distribution of C. jejuni CCs was different among continents, with ST-21 CC as the main one in North America, Europe, and Asia; ST-353 CC in South America and Africa, and ST-354 CC in Oceania. Moreover, this distribution by continent was not maintained among all the countries within a given continent, i.e., Denmark harbored a strong dominance of ST-21 CC genomes while Finland had ST-45 CC as the most abundant CC. Surprisingly, a lack of main association within developing-developed countries and AR composition on C. jejuni and C. coli genomes was clearly found, maybe due to the unequal distribution of both CC-STs and host origin of isolation among countries and even the amount of sequenced genomes for each country.
Results obtained in the current study reveal the association of certain C. jejuni CCs and both C. jejuni and C. coli STs to particular isolation sources. In this regard, it is worth mentioning the high abundance of CCs linked to chicken as a host. Most of these C. jejuni CCs had been already classified as a chicken specialist (Sheppard et al., 2014), and our study confirmed these previous associations, while ST-283 CC, previously described as a chicken specialist (Sheppard et al., 2014), can be considered as a Those CCs with less than 20 genomes with a known non-human source of isolation were not employed for the specialist/generalist host association analysis. The asterisk indicates CCs with new classification of host specialism/generalist behavior proposed, according to the data obtained in the current study. ST-22 CC is proposed to be cattle specialist instead of generalist and ST-283 CC to be generalist instead of chicken specialist (see Figure 4). to the results here reported. The high abundance of genomes from generalist ST-21 and chicken specialist ST-353 CCs agrees with previous reports that have described these CCs as highly prevalent among clinical, food, and animal samples from Israel (Rokney et al., 2018); clinical cases, dairy cattle, and broiler products in Lithuania (Ramonaite et al., 2017); chicken in Senegal (Kinana et al., 2006); or human and broiler carcasses in Belgium (Elhadidy et al., 2018). The high abundance of ST-21, ST-48, ST-45, and ST-257 CCs on isolates obtained from human samples has been previously highlighted in the United Kingdom (Haldenby et al., 2020). Interestingly, the STs within ST-45 CC, previously considered as a generalist CC, could be considered as chicken specialist STs. Conclusions on host specialization have to be taken carefully. For example, the generalist ST-21 CC harbored STs considered as chicken specialist (ST-50 and ST-53), cow specialist , and generalist . This host specialization differences within ST-21 CC could be due to the paraphyletic phylogeny of this CC, which was found as the ancestor CC of the cattle specialist ST-61 CC (Sheppard et al., 2014;Mourkas et al., 2020). The prevalence of ST-353, ST-354, ST-50, and ST-45 on chicken hosts and ST-21 as a host generalist were previously reported by using all the profiles on PubMLST , which are twice the number of genomes available on that web repository and employed in the current study. That supports the proposed host specialization classification, despite the possible biases introduced due to the isolates selected or employed for whole genome sequencing within all the available collections of bacterial isolates. Determinants of resistance to β-lactams and tetracyclines were the most extended among the 39,798 C. jejuni and 11,920 C. coli genomes analyzed, with bla OXA−61 as the main β-lactam resistance factor on C. jejuni; bla OXA−61 , bla OXA−489 , and bla OXA−193 on C. coli; and tet(O) for tetracycline resistance on both species. Determinants of resistance to quinolones [due mainly to the gyrA (T86I) SNP] and to aminoglycosides [by aph(3 )-III, aadE-Cc, and aph(3 )-VIIa genes, among others) were also highly prevalent on C. jejuni and C. coli genomes, respectively. Phenotypic resistance to β-lactam antibiotics in C. jejuni is frequently linked to the production of β-lactamase enzymes, encoded by bla OXA genes (Jonker and Picard, 2010), while tetracycline and quinolone resistance has been often associated with the tet(O) gene (Meistere et al., 2019) and the gyrA (T86I) SNP, respectively (Kinana et al., 2006;Meistere et al., 2019;Abraham et al., 2020;Joensen et al., 2020). Indeed, a high prevalence of bla OXA−61 and tet(O) has been found on C. jejuni isolated from human, chicken, and cattle in Spain in a study where SNPs were not screened (Mourkas et al., 2019) and on Campylobacter spp. in a global study analyzing 237 genomes from the RefSeq NCBI database (Rivera-Mendoza et al., 2020). Moreover, AR determinants of aminoglycoside resistance on C. coli and quinolone resistance on C. jejuni were mainly found on genomes associated with multidrug resistance, together with AR determinants for β-lactams and tetracyclines.
The increase in resistance to quinolones on Campylobacter isolates during the last years has been reported worldwide, including the United States (Nachamkin et al., 2002), China (Wang et al., 2016), France (Gallay et al., 2007), Vietnam (Nguyen et al., 2016), and Peru (Pollett et al., 2012), among other countries, and the usage of fluoroquinolones in poultry farms has been associated with an increase in resistance to quinolones in chicken and human Campylobacter isolates (Wieczorek and Osek, 2013). An increased occurrence of determinants of resistance to β-lactams, tetracyclines, and quinolones was found among those C. jejuni and C. coli genomes from human and food animal origin. The clear differences observed in the resistome between Campylobacter genomes from animals employed for food production and those from non-livestock animals can be likely related to the selective pressure exerted by antibiotic use in veterinary settings. Isolates from livestock species were characterized by a higher incidence of bla OXA−61 , associated with resistance to β-lactam antibiotics, such as amoxicillin or ampicillin, antibiotics that are frequently employed for veterinary uses. Quinolones and tetracyclines are still available for treatment of livestock all over the world (WHO, 2017), so this is likely the main cause of the resistance levels found for such antibiotics on isolates from food-related animals. Antibiotic resistance represents a great concern mainly in livestock animals due to the misuse of antibiotics. As an example, erythromycin (macrolide), ampicillin (β-lactam), tetracycline (tetracycline), nalidixic acid, and ciprofloxacin (both belonging to the quinolones family) are the drugs toward which the most resistant isolates were found in a recent review on Campylobacter spp. in Sub-Saharan African countries (Hlashwayo et al., 2020), which is in agreement with results found globally or also with the temporal changes among European isolates described in the current study.
Although the current study does not involve a source attribution analysis, the similar ARG genotype patterns and levels of prevalence of bla OXA−61 , tet(O), and gyrA (T86I) found among human, chicken, and turkey C. jejuni genomes, and of bla OXA−61 and tet(O) among human and chicken C. coli genomes, suggest a high level of relatedness among isolates from these origins. Poultry has been reported as the main source of human infection in studies carried out in different locations including Denmark (Joensen et al., 2020), Sub-Saharan Africa (Gahamanyi et al., 2020), or Lithuania (Ramonaite et al., 2017). On the other hand, other authors have identified as important sources of human infection both chicken and ruminants in France (Berthenet et al., 2019), the United States (Kelley et al., 2020), Israel (Rokney et al., 2018), and Sub-Saharan Africa (Hlashwayo et al., 2020) or chicken and wild birds in Baltic countries (Mäesaar et al., 2020). Moreover, the detection of no major differences in the resistance profiles among isolates from different points of the chicken meat processing chain (Dramé et al., 2020) suggests that AR spread is originated on upstream steps, such as the animal feeding, where the animals are exposed to antibiotics. A clear link between human and chicken Campylobacter isolates has been previously demonstrated by combining pulsed-field gel electrophoresis (PFGE) to subtype strains and antibiotic resistance phenotypic testing (Zhao et al., 2015).
The spread of AR genes among Campylobacter isolates from humans, animals, and the environment has been reported previously (Mourkas et al., 2019), and it should explain the presence of high levels of ARGs conferring resistance to tetracyclines in human isolates despite these antibiotics are not clinically employed in human medicine. Moreover, tetracycline is broadly employed in the poultry industry worldwide (Van Boeckel et al., 2015), which could be an important reservoir for resistant strains. Furthermore, higher levels of resistant C. jejuni FIGURE 5 | AR genes and SNP distribution among C. jejuni and C. coli genomes by source of isolation. Heatmap plot showing the percentage of (A) C. jejuni genomes and (B) C. coli genomes per isolation source carrying ARGs or SNPs associated with AR, showed in columns. Only the main 20 ARGs or SNPs are presented, because the clustering by isolation source remained the same as when the 120 and 68 AR determinants were used for C. jejuni and C. coli, respectively (data not shown). Only isolation sources containing more than 50 genomes are presented. Gene names are colored according to the corresponding antimicrobial family (blue: beta-lactams, orange: tetracyclines, red: aminoglycosides, purple: quinolones, and green: MLSP). Only the main 20 ARGs or SNPs are presented, because the clustering by isolation source remained the same as when the 120 AR determinants were used (data not shown). The Host_range classification of CCs is according to previous publications (see Table 1), except for those with highlighted edges, where edge color indicates previous classification while fill color indicates the new assignment, according to the results showed in Figure 4. Gene names are colored according to the corresponding antimicrobial family (blue: beta-lactams, orange: tetracyclines, red: aminoglycosides, purple: quinolones, and green: MLSP). and C. coli isolates were observed among isolates from pigs than from wild boars in Italy (Marotta et al., 2020), where an important amount of multidrug-resistant strains were isolated from pigs, being the over-use of antibiotics the main cause of the spread of multidrug resistance among pig isolates.
The bla OXA−61 gene is highly prevalent among Campylobacter spp. isolated from poultry (Griggs et al., 2009). The high diversity and abundance of determinants of resistance to β-lactams different to bla OXA−61 among isolates from nonfood-related animals and environments could be due to the ability of Campylobacter spp. for intrinsic production of β-lactamases in the absence of selective pressure (Cody et al., 2012). The bla OXA−61 gene, chromosomally encoded (Rivera-Mendoza et al., 2020), can generate many variants with few SNPs, differing from other bla oxa -encoded proteins in only one (e.g., bla OXA−193 ) or two amino acids (e.g., bla OXA−489 ). The bla OXA−61 gene and aminoglycoside resistance genes are confined to C. jejuni and C. coli genomes within Campylobacter species, while gyrA T86I SNP was mostly presented on C. jejuni (Rivera-Mendoza et al., 2020). Unfortunately, ACT to ATT mutation on gyrA position 86, the most common on C. coli strains with resistance to quinolones (Zirnstein et al., 2000), was not included on PointFinder database, which includes ACA to ATA mutation, the most prevalent on C. jejuni. So that, the prevalence of quinolone resistance among C. coli genomes found on the present study can be due to this technical issue.
The clonal complexes ST-354, ST-573, ST-464, and ST-446, the first three considered chicken specialist, were the CCs with the highest prevalence of ARGs linked to β-lactam, tetracycline, and quinolone resistance, due to the presence of the bla OXA−61 , tet(O), and gyrA (T86I) determinants for all of them. Some of these CCs  have been previously associated with ciprofloxacin (quinolone) resistance in human isolates (Cody et al., 2012). ARGs were more randomly distributed among CCs than among isolation sources. Hence, ARG genotype patterns were more closely related to the source of isolation (and probably the use of antibiotics) than to the lineage of the strains, as also remarked previously (Cody et al., 2012).
An increase in the occurrence of genotypes of resistance to tetracyclines and quinolones, mainly due to the spread of bla OXA−61 -gyrA (T86I)-tet(O) and bla OXA−61 -gyrA (T86I)tet(O/W/O) multi-drug resistance genotypes, was observed along the last 20 years in European C. jejuni genomes. The main reason for this temporal trend could be the continuous misuse of tetracyclines and quinolones, still broadly employed on livestock all over the world (WHO, 2017). Interestingly, fluoroquinolones are not utilized in Australian livestock, where Campylobacter isolates from chicken meat show very scarce levels of resistance to fluoroquinolones, always due to gyrA (T86I) SNPs (Abraham et al., 2020). On the other hand, an increase over time in the prevalence of the gyrA (T86I) SNP among isolates of C. jejuni associated with human gastrointestinal disease has been reported in the United Kingdom from 1993-1996-2009and even to 2016(Haldenby et al., 2020. This spread of resistance to quinolones is an important concern in some countries, such as Latvia, where resistance to quinolones was found, due to gyrA (T86I), in 90% of the C. jejuni isolates obtained from clinical and foodrelated sources from 2008 to 2016 (Meistere et al., 2019).
Remarkably, it has been demonstrated that high levels of fluoroquinolone resistance can persist in poultry even after discontinued use of these antibiotics (Price et al., 2005) and that resistance can rapidly emerge in poultry flocks (Humphrey et al., 2005). Some fitness benefits linked to the gyrA (T86I) SNP could be the reason for such high resistance prevalence in settings where the use of fluoroquinolones has been reduced, as has been proposed by Haldenby et al. (2020). Another example of temporal increases in fluoroquinolone resistance among Campylobacter species was observed on fecal samples from cattle in France between 2002 and 2006, with the fluoroquinolone resistance rate rising from 29.7% to 70.4% (Châtre et al., 2010).
In conclusion, the current work represents the biggest study to date, as far as we know, of the resistome markers harbored in C. jejuni/coli genomes, with 39,978 C. jejuni and 11,920 C. coli genomes in silico screened to search both for antimicrobial resistance genes and SNPs conferring antimicrobial resistance. The results obtained suggest the association between antimicrobial use in veterinary settings, in particular in poultry production, and the spread of resistance to Be, Qu, and Te and evidence the rapid expansion in Europe in the last two decades of determinants of resistance to these antibiotic families, especially among poultry and human isolates. The workflow employed is available to be used in future global resistome analyses for other species of interest, or even to be adapted to other kinds of microbiome analyses.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
AÁ-O and JC-D created the study design. JC-D and PG wrote and employed the Ruby and R scripts and downloaded genomes, made data analysis, and plotted the figures. JC-D and PG wrote the first draft of the manuscript and AÁ-O revised it. All authors revised the manuscript and approved the final version.

ACKNOWLEDGMENTS
We acknowledge the teams involved in the maintenance of PUBMLST, NCBI, and PATRIC web-sites, for making available such amount of information and data for the entire scientific community.