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¹; (ii) through the Export dataset section on the PubMLST-Campylobacter jejuni/coli website², where all isolates were selected; and (iii) by using 02.NCBI_Download_metadata.R script³ 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,⁴ PATRIC,⁵ and PubMLST⁶ 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⁷ 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).

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.

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.

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 ST-353 and ST-45 CCs (8.9 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).

The first global analysis of C. jejuni CC distribution showed differences by continent. The ST-21 CC was the most abundant 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).

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 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 (Table 2). According to the clustering obtained, ST-460, ST-464, ST-574, ST-607, and ST-658 CCs, previously not classified according to their host specialist/generalist behavior, can be classified as chicken specialist CCs, while ST-52 CC can be classified as a generalist CC (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 with medium levels of chicken specialization, while ST-353 and ST-354 CCs contain STs highly specialized for chicken (Figure 3B and Supplementary Figure 2B). Finally, the ST-21 CC contains three STs specialized in cattle infection (ST-8, ST-806, and ST-982), two specialized in chicken infection (ST-50 and ST-53), one ST specialized on sheep infection (ST-262), and two generalist STs (ST-19 and ST-21) (Figure 4B and Supplementary Figure 2B). The ST-858 CC, which was the main one found for C. coli genomes, harbors up to seven STs (ST-4425, ST-830, ST-1119, ST-872, ST-1614, ST-828, and ST-829) that could be considered as chicken infection specialist, three pig specialized STs (ST-1058, ST-854, and ST-887), one cattle specialist ST (ST-1068), and five host generalist STs (ST-1016, ST-1055, ST-827, ST-825, and ST-962) (Figure 4C and Supplementary Figure 2C). The complete list of STs analyzed, from both C. jejuni and C. coli genomes, and their host specilization were indicated on Supplementary 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).

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 clear association with the host range of each CC was observed at gene level.

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.8% in the 2000–2004 5-year period to 44.0% in the 2015–2018 period, and quinolones, from 21.0% in the 2000–2004 5-year period to 42.4% in the 2015–2018 period (Figure 7A and Supplementary Figure 5A).

These observations were linked to an increase in the occurrence of both the tet(O) and tet(O/W/O) genes and the gyrA (T86I) SNP (Figure 7B and Supplementary 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).