A total of 10 small specialty crop farms (SSCF) in Northeast Ohio (United States) were selected for this study, based on their availability for longitudinal sampling throughout the study period (2018–2020; Hailu et al., 2021). Farms were asked to continue performing management activities, which followed normal and customary practices of producers utilizing limited mechanization. The spatial distribution of the sample collection sites is provided in Figure 1. Each farm cluster (FC; n = 5) was created based on the geographic proximity of the farms with each other. The dairy manure was obtained from open dairy heifers housed on a bedded pack during the winter and raised on pasture in the summer. Dairy manure was not treated by the farmers. Dairy manure samples were collected monthly across the 10 SSCF between 2018 and 2020 (n = 140 total samples). Some samples were missing due to events surrounding the COVID-19 pandemic. Samples were collected aseptically into Nasco Whirl-Pak™ bags (Fisher Scientific, Waltham, MA, United States) and stored on ice until processed immediately in the lab.

Manure samples (25 g) were resuspended into 225 mL of phosphate buffered saline (PBS; Fisher Scientific, Waltham, MA, United States; 1:10 ratio), homogenized, serial diluted (10-fold), plated on RAPID’ Campylobacter (BioRad, Hercules, CA, United States), and then incubated at 42°C in microaerobic conditions (5% O₂, 10% CO₂, and 85% N₂) for 48 h. Resuspended manure samples (1 mL) were also enriched for 48 h at 42°C in 9 mL of Preston enrichment broth containing Campylobacter growth supplements (product No. CM067, Lysed Horse Blood product No. SR0048, Preston Campylobacter Selective Supplement product No. SR0117 and Campylobacter Growth Supplement product No. SR0232; Oxoid Ltd., Cambridge, United Kingdom; 1:10 ratio), plated on RAPID’ Campylobacter, and incubated under microaerobic condition at 42°C for 48 h. Colonies growing on RAPID’ Campylobacter plates were sub-cultured onto a fresh modified charcoal cefoperazone deoxycholate agar (mCCDA; Fisher Scientific, Waltham, MA, United States) plated and incubated at 42°C for 48 h under microaerobic conditions. A Campylobacter genus-specific colony-PCR was performed on all the Campylobacter-suspected colonies, as described below. Isolates confirmed to be Campylobacter were frozen at −80°C in Brucella broth supplemented with 30% glycerol (v/v).

Colony-PCR was performed on all isolates suspected to be Campylobacter using genus-specific primers (MD16S F: ATCTAATGGCTTAACCATTAAAC; C1228R: GGACGGTAACTAG TTTAGTATT; Final concentration = 0.11 μM; product size = 857 bp; Denis et al., 1999). Briefly, colonies were resuspended in 100 μL sterile DNase-free water and boiled at 95°C for 10 min. The product was centrifuged for 10 min at 4,000 ×g and 2 μL of the lysate was used for the colony-PCR (25 μL final volume; 25 cycles). The following program was used: one cycle of 10 min at 95°C, 35 cycles of 30s at 95°C, 1.5 min at 59°C, 1 min at 72°C, and a final extension step of 10 min at 72°C. PCR products were visualized using 2% agarose gels. Isolates confirmed as Campylobacter by PCR were selected for whole genome sequencing.

Genomic DNA from pure Campylobacter isolates (n = 69) was extracted using the Wizard Genomic DNA Purification kit (Promega, Madison, WI, United States). The DNA was quantified using a Qubit dsDNA BR assay kit (Invitrogen, Thermo Fisher Scientific, Waltham, MA, United States) and Qubit 3.0 fluorometer (Invitrogen). The extracted DNA samples were diluted in molecular biology grade nuclease-free deionized water (Thermo Fisher Scientific) to a final concentration of 0.2 ng/μL for whole genome sequencing (WGS). Genomic libraries were constructed using the Nextera XT DNA sample preparation kit as described in the manufacturer’s protocol (Illumina, San Diego, CA, United States) and sequenced on the MiSeq platform using 500 cycles of paired-end reads (Illumina). Quality control of the FASTQ datasets (raw reads) obtained from each sequence run was assessed using FastQC (Andrews, 2010). Low-quality reads < Q20 were trimmed and adaptor sequences were removed using the default parameter in Trimmomatic (Bolger et al., 2014). Trimmomatic modifications were conducted with a sliding window trimming option where the number of bases to average across was set to 4 and the average quality required was set to 20. The cleaned reads were subsequently de novo assembled using SKESA (Strategic K-mer Extension for Scrupulous Assemblies; Souvorov et al., 2018). Assembly qualities were also evaluated by QUAST (Quality Assessment Tool; Gurevich et al., 2013; Mikheenko et al., 2018). The analytical tools for FastQC, Trimmomatic, SKESA, and QUAST are available on the GalaxyTrakr platform (https://galaxytrakr.org, accessed 27 May 2022) which is an open-source bioinformatics platform maintained by the Center for Food Safety and Applied Nutrition (CFSAN) at the U.S. Food and Drug Administration (Gangiredla et al., 2021).

Sequence types (STs) of Campylobacter genomes investigated in this study were initially analyzed using the traditional seven-loci Multi-locus sequence typing (MLST) scheme based on the following housekeeping genes (aspA, glnA, gltA, glyA, pgm, tkt, and uncA; Dingle et al., 2001). Briefly, the genome assemblies were scanned against the PubMLST schemes where the ST of each genome was assigned by comparing the alleles of the seven genes in the MLST open database (Jolley and Maiden, 2010; Seemann, 2022b).

To cover the expanded resolution of the MLST concept, a core genome MLST (cgMLST) profiling against the PubMLST database (Jolley et al., 2018) was performed using a modified Seemann’s MLST tool which used contig files to scan against the Oxford C. jejuni and C. coli cgMLST scheme (Mulder et al., 2020). The cgMLST scheme describes the genetic variation among the strains by utilizing the allele sequences of 1,343 loci which were defined from genome sequences of 2,472 representative United Kingdom campylobacteriosis isolates including C. jejuni and C. coli (Cody et al., 2017). The traditional MLST tool was used on the GalaxyTrakr platform (accessed 28 April 2022) while the cgMLST tool was accessed (28 April and 16 December 2022) from an open-source platform¹ to query the genomes against the current C. jejuni and C. coli scheme. The closest ST match for these isolates were assigned using the cgMLST scheme on the PubMLST Campylobacter jejuni/coli typing database² (Cody et al., 2017; Jolley et al., 2018).

To identify the isolates at the genus and species level, the sequence data were analyzed by Kraken 2, which is a k-mer-based taxonomic sequence classifier using k-mers of 35 bp (Wood and Salzberg, 2014; Wood et al., 2019). The Kraken 2 algorithm assigns taxonomic labels by matching each k-mer within a query sequence to the lowest common ancestor (LCA) of the pre-computed database of genomes containing the given k-mer (Wood et al., 2019). In addition, the genomes were annotated by uploading the assemblies of FASTA datasets through the National Center for Biotechnology Information (NCBI) prokaryotic genome annotation pipeline (PGAP) with its best-placed reference protein set GeneMarkS+ application (Haft et al., 2018; Li et al., 2021). Assembled genomes were also annotated using RAST (Aziz et al., 2008; Overbeek et al., 2014) to provide an overall picture of the gene content across the Campylobacter isolates.

The predicted antimicrobial resistance (AMR) genes were identified using the National Antimicrobial Resistance Monitoring System (NARMS) Campylobacter workflow on GalaxyTrakr which uses BLAST techniques against the NCBI’s comprehensive AMR gene database (Feldgarden et al., 2019a,b, 2021). The AMRFinderPlus tool accessed from the GalaxyTrakr platform was used to detect acquired AMR genes in bacterial proteins or assembled nucleotide sequences, along with point mutations that were cross-referenced with a core set of AMR elements and the expanded subset of reference genes related to biocide, stress response, and virulence factors (Feldgarden et al., 2021). The AMRFinderPlus tool was designed to utilize the Reference Gene Catalog database of NCBI Pathogen Detection, which consists of 6,428 genes (5,588 AMR genes, 210 stress response genes, and 630 virulence genes), 627 hidden Mark models (HMMs), and 682 point-mutations (Feldgarden et al., 2021). Virulence gene profiles of Campylobacter isolates were determined by BLASTN comparison against the virulence factor database (VFDB) with the threshold of 80% identity and 80% sequence coverage (Chen et al., 2016; Liu et al., 2019) by using an ABRicate tool which performed mass screening of contigs for antimicrobial resistance and virulence genes (Seemann, 2022a).

To understand relatedness among the strains, the Campylobacter genome assemblies were analyzed using the CFSAN SNP Pipeline (Davis et al., 2015) which is accessible on the GalaxyTrakr platform. The single nucleotide polymorphism (SNP) pipeline performs the alignment of mapped reads of sequences to a reference genome, which creates high-quality SNP matrices for sequences and determines phylogenetic relatedness among the strains (Davis et al., 2015). The SNP analysis was conducted against the reference genome chosen by the SNP pipeline, and a reference strain was chosen based on the best assembly metrics from the quast report, such as the one with the longest N50. For this analysis, sample #112962 (Biosample accession # SAMN29955749) was selected as a reference to build the tree. The phylogenetic tree was constructed using the Neighbor-Joining method (Saitou and Nei, 1987). The evolutionary analyses were employed using MEGA7 software (Kumar et al., 2016, p. 7), with evolutionary distances computed using the Maximum Composite Likelihood method (Tamura et al., 2004).

Assembled genomes were scanned for plasmids using PlasmidFinder v. 2.1.6 (Carattoli et al., 2014), which detects plasmids via replicons only, and MOB-Suite v. 3.1.0 (Robertson and Nash, 2018), which detects plasmids with a database containing replicons, relaxases, and complete plasmid sequences. PlasmidFinder was run using the plasmidfinder.py script, after downloading the database with the download-db.sh script (on 16 July 2022), with a minimum coverage threshold of 0.60 (the default) and an identity threshold of 0.95 (the default for the webserver, whereas the default for the Python script is 0.90). MOB-Suite was run via the “mob_recon” command, after running “mob_init” to download the databases (on 18 November 2022) with default parameters.

Genomic data were synchronized with the metadata (e.g., collection time point and location, and source) and analyzed using JMP Pro 16 software (SAS Institute Inc., Cary, NC, United States). The analyses described below focus essentially on the variable genomes (protein-encoded genes not consistently detected within all the C. jejuni or C. coli genomes) and not the core genome (protein-encoded genes consistently detected within all the C. jejuni or C. coli genomes). Multi-dimensional scaling analysis combined with K-means clustering were used to create a two-dimensional plot displaying the farm geographic distribution and associated clusters based on the GPS coordinates of the SSCF. Hierarchical clustering analysis was conducted to identify similar variable gene profiles among the 69 C. jejuni and C. coli isolates. Principal component (PCA) and agglomerative hierarchical clustering (HCA) analyses were used to study the distribution of the Campylobacter isolates based on their functional genomic profiles and cgMLST. Multi-correspondence analysis (MCA) combined with a Chi² test were used to determine whether the prevalence of certain genes was influenced by the Campylobacter species, SSCF, and sample collection time point (e.g., year and month of collection). MCA and multivariate analyses were used to identify co-occurrence of specific genes for a given Campylobacter species. In addition, an association study (also called produce market analysis) was conducted to determine whether the prevalence of specific genes could be linked with specific geographic or temporal parameters associated with the isolates. A minimal confidence level of 95% and lift of 2 were used to identify associations of statistical importance.

Out of the 140 dairy manure samples collected from the 10 SSCF between 2018 and 2020, 49% (n = 69) were positive for thermophilic Campylobacter (Table 1). Dairy manure collected in 2020 and 2019 harbored the highest prevalence (88% and 69%, respectively) compared to 2018 (18%; p < 0.001). On the other hand, dairy manure collected in 2020 harbored the lowest abundance of Campylobacter (2.92 [IC95%: 2.61–3.22]) log CFU/g compared to 2019 and 2018 (3.89 [IC95%: 3.64–4.13] and 3.28 [IC95%: 2.99–3.57] log CFU/g, respectively Table 1; p < 0.008). Equivalent Campylobacter prevalence (60%–100%) and load (2.76 [IC95%: 2.31–3.21] to 4.21 [IC95%: 3.82–4.60] log CFU/g) in the dairy manure were detected between farms and seasons (p > 0.01).

The whole genome analyses revealed 56 C. jejuni (genome length = 1,683,112 bp [IC95%: 1,664,247–1,701,977]) and 13 C. coli (genome length = 1,691,610 bp [IC95%: 1,680,103-1,703,117]) isolates (Supplementary Table 1). All associated metadata and genome quality statistics are included in Supplementary Table 2. Campylobacter jejuni was recovered across all 3 years of sampling, predominantly in 2019 (n = 33/56; Table 2). Most of the C. coli isolates were recovered in 2020 (n = 12/13; Table 2). Isolates were collected in the winter, spring, and summer seasons (n = 16–19 C. jejuni and 3–6 C. coli isolates per season), but rarely in the fall (n = 3/56 C. jejuni isolates and no C. coli isolates).

A total of 20 STs were identified among the 56 C. jejuni isolates, with ST-922 (n = 10) and ST-61 (n = 7) being the most predominant (Table 3). These 20 STs belonged to 10 clonal complexes (CCs), with CC-21 (n = 26) being the most predominant. The geographic location of the farms affected the type of ST and CC detected. The majority of the ST-61 (n = 6/7) were detected in farm #30 in 2019; ST-922 (n = 7/10) were detected in farm #26 in 2019 and ST-829 (n = 5/8) were detected in farm #7 in 2020. The spatiotemporal distribution of the CC and ST was studied across the 10 farms over 3 years. Similar CCs were detected in FC #1 (farm #1, 7 and 26) and FC #2 (farm #8 and 1 7; Supplementary Figure 1A), which is in concordance with the geographic proximity of both FC (Figure 1). However, FC #3, #4, and #5 harbored different CCs compared to FC #1 and #2 (Supplementary Figure 1A), which is also in agreement with the distant geographic location of the FCs between each other (Figure 1). Two STs (ST-829 [n = 8] and ST-1068 [n = 5]) were detected among the 13 C. coli isolates and belonged to CC-828. The number of STs and CCs (C. jejuni and C. coli combined) were equivalent across collection years (n = 9–12 and 5–7, respectively), but each year a different ST/CC profile was observed (Supplementary Figure 1B). No differences in ST and CC profiles were detected between seasons.

Distinct SNP signatures were detected between different Campylobacter species, as well as, between isolates of the same species (Figure 2). As expected, C. jejuni isolates segregated from the C. coli isolates. A total of six major C. jejuni clusters were observed based on the SNP profiles. The composition of these clusters was associated to some extent with the geographic location, and less with the time of collection (Figure 2). Most of the C. jejuni isolates collected from farm #30 (n = 6/8) clustered together and displayed identical SNP profiles (cluster 2). The majority of C. jejuni isolates collected from FC #1 (n = 16/21; farm #1, #7, and #26; Figure 1) clustered away from the other C. jejuni isolates (n = 31; cluster 4; Figure 2). Interestingly, FC #1 was mainly composed of C. jejuni isolated between the summers of 2018 and 2019, and harbored high SNP profile similarity with a C. jejuni isolated in early 2018 from farm #11. Farm #11 belongs to FC #4 which is in close proximity to FC #1 (Figure 1). High SNP profile similarities were also observed between the C. jejuni isolates collected between FC #1 and FC #2 (farm #8 and #17; cluster 3 and 5; Figure 2). A similar trend was observed with C. coli isolates between FC#1, FC#2 and FC#5 (C. coli cluster; Figure 2).

The 69 Campylobacter isolates were also analyzed using core genome MLST (cgMLST). A total of 869/1,343 alleles were detected across the 69 isolates (Supplementary Figure 2). Overall, a distinct separation between C. jejuni and C. coli was observed. For the C. coli isolates, ST-829 isolates (n = 8) clustered away from the ST-1068 isolates (n = 5). Interestingly, C. jejuni CC-21 isolates were subdivided into four clusters; all C. jejuni ST-922 CC-21 isolates (n = 10) clustered away from the second (i.e., C. jejuni CC-21 ST-8 [n = 5]), third (i.e., C. jejuni CC-21 ST-929 [n = 3]) and fourth cluster (i.e., C. jejuni CC-21 ST-806 [n = 4] and ST-21 [n = 1]). The allelic profile of C. jejuni ST-922 CC-21 isolates was also closely similar to C. jejuni CC-403, CC-463 and CC-52. Similarly, other C. jejuni CC-21 (i.e., ST-806 and ST-982) were closely related to C. jejuni CC-48 and CC-61. Campylobacter jejuni ST-5261 CC-257 were closely related to of C. jejuni ST-929 CC-21 cluster, which itself displayed similarities with C. jejuni CC-42 isolates. Only C. jejuni ST-8 CC-21 and ST-45 CC-45 isolates stood apart from other C. jejuni isolates.

After annotating the genomes using RAST (general gene content profile), VFDB (virulome annotation) and NDARO (resistome annotation), a total of 1,784 annotated genes were detected across the 69 Campylobacter isolates. Among them, 35.7% (n = 637/1,784) of annotated genes were not always detected in the 69 Campylobacter isolates studied (referred to as variable genes). As observed with the SNP analysis (Figure 2), C. coli isolates displayed a different gene content profile compared to the C. jejuni isolates. Campylobacter jejuni isolates possessed several unique virulence-associated genes involved in the uptake of potassium, carbon utilization (e.g., succinate, gluconate, oxoglutarate and malate), and chemotaxis (n = 35; Supplementary Table 3). Similarly, C. coli isolates also had unique virulence genes associated with capsule production, cell envelope integrity, and membrane transporters, increasing resistance to environmental stresses (Supplementary Table 3). Campylobacter coli isolates also possessed unique sets of genes coding for enzymes (e.g., histidine kinase, acylamide amidohydrolase, methylcitrate dehydratase and synthase, adenine-specific DNA methyltransferase, adenylylsulfate kinase, L-carnitine dehydratase, and methionine synthase), transporters (e.g., creD, yihN, and yihY) essential for the utilization of specific energy sources, and proteins involved in iron availability and/or acquisition (e.g., ycsG, hemerythrin, hemolysin, and peroxide stress regulator/ferric uptake regulation protein; Supplementary Table 3).

A total of 401 protein-encoding genes were identified as part of the C. jejuni variable genome (all variable genes identified across the 56 C. jejuni genomes). Specific gene content profiles (n = 12 major clusters; Table 4; Supplementary Table 4) were identified. We observed that the gene content profiles were closely associated with ST of the isolates (Figure 3). Additional details about the variable gene profile distribution observed for each C. jejuni isolate are presented in Supplementary Figure 3A.

Most C. jejuni isolates (n = 51/56 isolates) possessed genes associated with lipooligosaccharide synthesis, O-methyl phosphoramidate capsule modification and methionine metabolism (n = 12 genes; Cluster Cj1 in Table 4; Supplementary Table 4). On the other hand, five isolates missing the genes from cluster Cj1 possessed genes involved in fatty acid synthesis, dmsABC sulfoxide reductase and capsular polysaccharide biosynthesis (n = 8 genes; Cluster Cj2 in Table 4; Supplementary Table 4). Additionally, most isolates (n = 52/56 isolates) possessed genes associated with pantothenate and rhamnose synthesis, capsular polysaccharide biosynthesis, fucose utilization, and bacteriocin resistance (n = 12; Cluster Cj3 in Table 4; Supplementary Table 4). However, the four isolates missing the genes from cluster Cj3 possessed other genes involved in arsenate resistance and glycosaminoglycan biosynthesis (n = 5 genes; Cluster Cj4 in Table 4; Supplementary Table 4).

A fraction of the isolates (n = 38/56 isolates) possessed genes associated with cell surface glycoconjugates (glutamate and legionaminic acid; n = 8 genes; Cluster Cj5 in Table 4; Supplementary Table 4). Among these isolates, 22 C. jejuni (especially from FC #2) also harbored genes for the utilization of fucose and glucuronate synthesis (n = 6 genes; Cluster Cj6 in Table 4; Supplementary Table 4).

Several isolates from FC #1 (n = 10/22 isolates; especially ST-922 CC-21 from farm #26) harbored genes associated with rhamnose, methionine, and purine metabolism (n = 10 genes; Cluster Cj7 in Table 4; Supplementary Table 4). The majority of these isolates also possessed type IV secretion systems (T4SS) protein-encoding genes (virB2,3,5–9, and virD4; n = 10 genes), IncQ plasmid conjugative transfer-related genes (traC, traE, traG and traR; n = 4 genes), and the gene encoding cag pathogenicity island protein (cag12; Cluster Cj8 in Table 4; Supplementary Table 4). A similar T4SS/conjugation profile was observed in C. jejuni isolates from FC #2 (n = 6/17 isolates, especially farm #17). These isolates also harbored genes associated with IncF plasmid conjugative transfer pilus assembly (n = 11 genes; traB,C,E,H,K,L,N,T,U,V, and W; Cluster Cj9 in Table 4; Supplementary Table 4). Further, 5 of these isolates also possessed additional protein-encoding genes associated with T4SS (virB4,7) and conjugative transfer (traG), hypothetical pVir protein (n = 10 genes, pVir0004,7,8,9,12,15,19,20,29,42), and other virulence-associated genes (n = 15; Cluster Cj10 in Table 4; Supplementary Table 4).

Another group of isolates (n = 12 isolates) harbored genes associated with conjugative transfer (n = 6; trbB,D,E,F,I,L) and T4SS (virB1), as well as a high abundance of bacteriophage DNA (n = 29 genes) especially from the Escherichia coli phage protein families Gp (phage lambda) and Mup (phage Mu; Cluster Cj11 in Table 4; Supplementary Table 4). Three isolates harbored type VI secretion system genes (T6SS; n = 6 genes; hcp, impB,C,G,I,K; Cluster Cj12 in Table 4; Supplementary Table 4).

Eighty-one protein encoding genes belonging to the variable genome were not clustered based on the CC or ST of the C. jejuni isolates. A majority of the protein-encoding genes (n = 27) were labeled as possible, putative, uncharacterized or hypothetical proteins (Supplementary Table 4).

A total of 102 protein-encoding genes were identified as part of the C. coli variable genome. Specific gene content profiles (n = 2 major clusters) were identified based on the ST associated with the isolates (Table 5; Supplementary Table 5). Additional details about the variable gene profile distribution observed for each C. coli isolate are presented in Supplementary Figure 3B. Interestingly, the major difference between these two clusters is the presence/absence of genes involved in the release and uptake of DNA and protein. All C. coli ST-829 CC-828 (n = 8 isolates) possessed phage DNA (n = 24 genes, especially from Gp and Mup families) and genes encoding enzymes associated with hexose utilization (glucose, abequose, fucose, mannose, glycerate, glucoronate and rhamnose; n = 14 genes; cluster Cc1 in Table 5; Supplementary Table 5). On the other hand, C. coli ST-1068 CC-828 (n = 5 isolates) possessed genes associated with conjugation transfer (n = 4 genes; traC,E,G, and Q), T4SS (n = 10 genes; virB2,4–11 and virD4), capsule modification (n = 5 genes), and cell wall synthesis (n = 2 genes; cluster Cc2 in Table 5; Supplementary Table 5).

Using the National Antimicrobial Resistance Monitoring System for Enteric Bacteria (NARMS) database, a total of 11 antimicrobial resistance genes (ARGs) were detected across the 69 Campylobacter isolates (Table 6; Figure 4). Three ARGs (blaOXA-193, gyrA_T86I, and tetO) were detected in both C. jejuni and C. coli isolates (Table 6). Campylobacter jejuni isolates were characterized with beta-lactamase (blaOXA-449, blaOXA-461, blaOXA-603, and blaOXA-61), macrolide (50S_L22_A103V), and aminoglycoside ARGs (rpsL_K88R and aph3’IIIa; Figure 4A), while C. coli isolates were characterized with an aminoglycoside ARG (aadE-Cc; Figure 4B). Interestingly, all C. coli isolates also harbored genes encoding proteins associated with resistance to arsenic (arsenate reductase [EC 1.20.4.4] thioredoxin-coupled and arsenical-resistance protein ACR3). All isolates with aadE-Cc also harbored tetO and gyrA_T86I (r² = −0.85; p = 0.0002), and all isolates with tetO also had gyrA_T86I (r² = 1; p < 0.0001) in C. jejuni. No ARG co-occurrence was detected in C. coli. A total of 14 ARG profiles were observed across the 69 isolates (Supplementary Table 6). Most of the isolates possessed genes conferring resistance to at least two antibiotics (n = 43/69 isolates). The co-occurrence of genes involved in the resistance to streptomycin (aadE-Cc) and beta-lactam (blaOXA-193) or to tetracycline (tetO), quinolone (gyrA_T86I), and beta-lactam (blaOXA-193) were predominant in C. coli (n = 7/13 and 5/13 isolates, respectively; Supplementary Table 7). Similarly, for most C. jejuni, the resistance to tetracycline (tetO) and beta-lactam (blaOXA-193, −461 or −449; n = 16/56 isolates), or to kanamycin (aph[3’IIIa]), beta-lactam (blaOXA-193), and tetracycline (tetO) were predominant (n = 8/56 isolates; Supplementary Table 7). Interestingly, gyrA_T86I and tetO were more likely detected in the C. coli isolates collected in the winter (60% confidence level and lift of 2), while aadE-CC was more likely detected in the C. coli isolates collected in the summer (57% confidence level and lift of 2.1). A similar trend was detected with blaOXA-193 in the summer with C. jejuni (40% confidence level and lift of 1.2).

The co-occurrence of genes from the variable genome was investigated in both C. jejuni (n = 401 genes) and C. coli (n = 102 genes). Overall, the C. coli variable genome was divided into three major clusters with co-occurring genes (Figure 5A; p < 0.01). The first co-occurring cluster named “conjugation/T4SS” (r² > 0.8; n = 47 genes; Supplementary Table 8) was composed of genes involved in conjugative transfer and T4SS (n = 17 genes), capsule associated genes (n = 5), cag12 (encoding cag pathogenicity island protein) and other virulence genes (e.g., rfbA, rfbB, and rhmA). Campylobacter coli isolates possessing the tetO gene, conferring tetracycline resistance, also possessed genes associated with conjugation/T4SS (r² = 1; p < 0.001). A similar observation was made regarding C. jejuni (r² > 0.43; p < 0.01). Interestingly, C. coli isolates possessing these conjugation/T4SS associated genes were less likely to possess prophages (n = 22 genes), and virulence genes (e.g., ceuB, fabG, glxK, rfbC, and rfbF; r² < −0.91; p < 0.00 l n = 47 genes; Figure 5A; Supplementary Table 8). The opposite trend between the co-occurrence of prophage and conjugative/T4SS clusters was observed in C. jejuni (r² = 0.69; p < 0.01). Both C. coli and C. jejuni isolates possessing genes implicated in O-methyl phosphoramindate capsule modification (hddA, hddC, and hddD) were also more likely to possess prophages (r² = 0.78 in C. coli and r² = 0.37 in C. jejuni; p < 0.01). Campylobacter coli isolates possessing genes associated with arsenical resistance (n = 3 genes) were likely to have genes involved in streptomycin resistance (n = 1 gene) and other genes (e.g., yrrC, yraQ family and ydeQ/yrkL/ywrO family). Unlike C. coli, the C. jejuni variable genome displayed less pronounced and indistinct gene co-occurrence profiles (Figure 5B).

While no plasmids were detected by PlasmidFinder, MOB-Suite identified a total of 223 plasmids across 69 assembled Campylobacter genomes. Only five plasmids were detected in C. coli (on average 0.54 plasmids per assembly), whereas 216 plasmids were detected in C. jejuni (on average 3.86 plasmids per assembly). Multiple plasmids were detected in each farm. Most plasmids were highly similar and grouped into a total of eight clusters (an approximation of Operational Taxonomic Units for plasmids) and only 18 unique reference plasmids were the nearest neighbors to each detected plasmid (Supplementary Figure 4; Supplementary Table 9). However, it is important to mention that sequencing was conducted using a protocol designed to extract genomic DNA and not plasmid DNA. By consequence, results presented in this section might underestimate the plasmid composition in the Campylobacter isolates. Overall, among the 18 plasmids, 10 of them harbored ARG associated with the resistance to tetracyclines and aminoglycosides, four of them harbored genes associated with the T6SS and two other plasmids harbored genes associated with the T4SS. The prevalence of certain plasmids was closely associated with the ST-CC profile of the isolates.