All Campylobacter jejuni and C. coli isolates used in this study were collected between December 2012 and December 2013 as part of the National Microbiological Baseline Study (MBS) conducted by the Canadian Food Inspection Agency in collaboration with industry and government partners (Canadian Food Inspection Agency, 2016). All isolates used in the study were recovered from both federally registered chicken establishments and retail outlets. Antimicrobial susceptibility testing (AST) was conducted as described in below. A subset of 271 isolates (229 C. jejuni and 42 C. coli) was selected for this study (Supplementary Table 1). The panel consisted of 109 isolates that were phenotypically sensitive to azithromycin (AZM), ciprofloxacin (CIP), erythromycin (ERY), gentamicin (GEN), tetracycline (TET), florfenicol (FLR), nalidixic acid (NAL), telithromycin (TEL), and clindamycin (CLI) and 162 isolates that were resistant to at least one of these antibiotics. Although none of the isolates used in this panel were phenotypically resistant to telithromycin, florfenicol, or gentamicin, the panel did include 37 isolates with an intermediate level of resistance to telithromycin (Supplementary Table 1).

In the case of disagreements between the AST profile and the WGS-derived profile, the AST profile of the isolate was verified as described in Dramé et al. (2020). Briefly, the minimum inhibitory concentration (MIC) for each of the nine antibiotics was determined by the broth dilution method in a microtitre plate format (Public Health Agency of Canada, 2015) and antibiotic susceptibility breakpoints were as described in Dramé et al. (2020) (Supplementary Table 1).

The AST was verified for 13 isolates by broth dilution, based on the methods described by the Public Health Agency of Canada (2015) and National Antimicrobial Resistance Monitoring System (2016) (Supplementary Table 2). Briefly, isolates were streaked for purity on Brucella agar from frozen stocks stored in Brucella broth (Oxoid) with 20% glycerol at −80°C. Isolates were cultured on Mueller Hinton agar (BD Difco) with 5% sheep blood (Cedarlane) at 42°C for 24 h in a microaerobic atmosphere. Cell were harvested from agar plates and cell density was to adjusted to a 0.5 McFarland standard in cation-adjusted Mueller Hinton broth with TES buffer (CAMHB+TES) (Thermo Scientific) using a DensiCHEK Plus instrument (bioMérieux). The cell suspension was then diluted 1:100 in CAMHB+TES with sheep blood (Thermo Scientific) and used as the inoculum for a broth dilution assay. The range of antibiotic concentration tested were as described by Public Health Agency of Canada (2015). Cultures were incubated in the microaerobic atmosphere at 42°C for 24 h. The minimum inhibitory concentration (MIC) was the lowest concentration at which no growth occurred, break-points for each antimicrobial were as described by Public Health Agency of Canada (2015) (Supplementary Table 1).

Isolates were grown from frozen glycerol stocks by streaking for single colonies on 5% horse blood agar. Cultures were incubated for 24 to 48 h at 42°C in microaerobic conditions. DNA was extracted from single colonies using the EZ1 DNA tissue kit (Qiagen) according to the manufacturer’s instructions. Genomic DNA was quantified using the Quant-it High-Sensitivity DNA Assay Kit (Life Technologies Inc., Burlington, ON, United States) and DNA libraries were prepared from 1 ng of genomic DNA with Illumina Nextera XT DNA sample preparation kit (Illumina, Inc., San Diego, CA, United States) according to manufacturers’ instructions. Whole-genome sequencing was performed on the Illumina MiSeq platform (Illumina Inc.), using a 600 cycle MiSeq reagent kit (v3) generating 2 × 300 bp paired-end reads.

Basic metrics of FASTQ-formatted sequencing reads were determined with FastQC version 0.11.6 (Andrews, 2010). Contamination was detected in sequencing reads using ConFindr (Low et al., 2019). Quality and adapter trimming of reads was performed with the BBduk script from the BBTools suite version 37.78 (Bushnell, 2014) with the following parameters: trim quality of 20, and remove reads below 50 bp long. Error correction was performed using tadpole version 37.78 (Bushnell, 2014) in “correct” mode with default parameters. Trimmed and corrected reads were assembled using SKESA v2.3.0 (Souvorov et al., 2018) and SPAdes v3.13.0 (Bankevich et al., 2012) with shovill 1.0.4¹. Pilon version 1.2.2 (Walker et al., 2014) was used to automatically improve assembly, and assembly metrics were calculated with a custom Python script², and Qualimap version 2.2.2 (Okonechnikov et al., 2016). Contigs shorter than 500 bp were removed from assemblies. Sequencing was repeated for strains with poor assembly metrics (e.g., low N50, coverage) or when contamination was detected.

Raw data have been deposited at DDBJ/EMBL/GenBank under BioProject PRJNA600010. The accession numbers and strain descriptions are listed in the Supplementary Table 1.

Identification of antibiotic resistance genes (ARGs) in the raw reads and draft genome assemblies (contig files) of each isolate was completed using the bioinformatics tools and databases listed in Table 1. Genes were deemed to be present based on a minimum sequence identity of 90%. Currently, PointFinder is unable to detect the gyrA gene of C. coli and at the time analysis the program was unable to reliably detect mutations in the 23S rRNA gene. As such, known point mutations in these two genes were screened for manually. Point mutations conferring AMR in the gyrA and the 23S rRNA genes were identified as follows. The gyrA and 23S rRNA gene sequences were obtained from the assemblies using GeneSeekr³ or from the raw reads using Sipprverse (Lambert et al., 2015; Table 1). Reference sequences for the target genes were from C. jejuni NCTC 11168 (accession number AL111168) and C. coli RM 4661 (accession number CP007181). Query gene sequences were identified by GeneSeekr through BLASTn (Altschul et al., 1990) sequence comparisons against the reference genes, and enabling the “–fasta_output” option. Sipprverse mapping of raw FASTQ reads to reference gene sequences with a kmer size of 11 and soft-clips enabled, allowed the extraction of the gene sequences. Sipprverse read-baiting was performed with BBduk version 37.78 (Bushnell, 2014). FASTQ reads were mapped to reference files using BOWTIE2 version 2.3.5 (Langmead and Salzberg, 2012), with “–local” alignment, “–very-sensitive-local” scoring, and “–all” alignments reported. BAM files were indexed and sorted with samtools version 1.9 (Li et al., 2009). Alignments were parsed from sorted BAM files using pysam version 0.15.2⁴. Sequence alignments of reference and query genes were done using Geneious Prime (version 2019.2) with nucleotide position numbering of the 23S rRNA gene based on the C. jejuni NCTC 11168 reference sequence. Mutations in the gyrA gene (amino acid position 86, Thr to Ile) which confer resistance to quinolone and the 23S rRNA gene (positions 2074, A to T, and 2075, A to G), which confer resistance to macrolides/lincosamides/ketolides, were identified.

For isolates with an identified 23S rRNA gene mutation (n = 42), sequences of cmeR, cmeABC operon, L4 (rspD), and L22 (rspV) were investigated for potential mutations relative to reference sequences from C. jejuni and C. coli. In addition, the presence of base substitutions (relative to C. jejuni NCTC 11168) in the CmeR binding site identified by Lin et al. (2005) were also investigated. Sequences for the genes of interest, including the tet(O) gene for isolates carrying the gene, were obtained using Sipprverse and aligned against reference sequences using Geneious (as described above). Isolates carrying a 23S rRNA gene mutation were screened for additional AMR markers using the Comprehensive Antibiotic Resistance Database (CARD) (Alcock et al., 2020).

To evaluate the accuracy of using WGS data to predict AMR in Campylobacter the following calculations were made using the phenotype profile as the reference method; overall association ([(a+d)/(a+b+c+d)]x100)), sensitivity ([a/(a+c)] × 100), specificity ([d/(b+d)] × 100), positive predictive value (PPV) ([a/(a+b)] × 100), and negative predictive value (NPV) ([d/(c+d)] × 100), where a = true positives, b = false positives, c = false negatives, and d = true negatives.

An estimate of the average genome coverage needed to reliably detect an ARG was determined by screening at increasing levels of read coverage. For each isolate, the raw reads were randomly sampled to genome coverage levels of 1, 2.5, 5, 10, 15, 20, 25, and 30X, using the reformat.sh script (version 37.61) provided with the BBMap suite (Bushnell, 2014). Each set of subsampled raw reads was then screened for ARGs using sipprverse (v 0.0.76) with the ResFinder database (2019-04-26), and SRST2 (v 0.2.0) (Inouye et al., 2014), and KMA (v 1.3.3) (Clausen et al., 2018) with the NCBI-AMR database (2019-04-29.1) (Table 1). For all three tools, default settings were used, with a minimum gene coverage was set to 90%. For each isolate, at each genome coverage level, the subsampling and ARG screening was repeated 100X. The expected ARG profile of each isolate was determined by running each of the tools against the raw reads without subsampling. A false positive result was defined as an ARG detected in the subsampled raw reads but not detected using all available raw reads.

Overall, there was a high concordance between the phenotypic susceptibility and the corresponding genotype (96.3–100%) (Table 2). Of the 271 isolates tested, discrepancies were observed in 13 isolates which were retested using the broth dilution method (Supplementary Table 2). Most of the discrepancies were due to discordance between the CLI and/or TEL genotype and the reported phenotype, i.e., presence or absence of the A2075G 23S rRNA gene mutation but lack of the expected susceptibility profile (Supplementary Table 2). Repeat testing led to the amendment of AST profiles for three of the isolates. Isolates CJ-MBS3796A and CC-MBS0767R were amended as being susceptible to CLI and CJ-MBS3139A was amended to being susceptible to AZM, ERY, TEL, and CLI (Supplementary Table 2). The amended phenotypes were used in subsequent analyses. The remaining isolates which were phenotypically susceptible to CLI and/or TEL but carried the A2075G 23S rRNA gene mutation, were confirmed as being susceptible to the antibiotic.

There was a 99.6% association between tetracycline resistance and the detection of a tet(O) determinant. The only discordance was a single isolate (CC-MBS7487A), in which a tet(O) gene was detected by all of the pipelines but the isolate was susceptible to TET (Supplementary Table 2). The sequences of the tet(O) gene from this isolate and all other isolates positive for a tet(O) gene (n = 128) were aligned against a tet(O) reference sequence (accession number M18896). The tet(O) sequence of CC-MBS7487A was the only sequence to have a substitution at nucleotide position 1609 (G → T) which resulted in a premature stop codon at amino acid residue 537.

Point mutation in the gyrA (T86I) or the 23S rRNA (A2075G) gene was found to correspond 100% to phenotypic resistance to quinolones (CIP and NAL) and macrolides (AZI and ERY), respectively (Table 2). In contrast, there was a lower association between a A2075G 23S rRNA gene mutation and increased resistance to TEL and CLI, with an overall 96.3 and 98.9% association between the predicted AMR and its phenotype for TEL and CLI, respectively. For both TEL and CLI, the ability of WGS to accurately predict sensitivity to each antibiotic was 100%, (i.e., isolates that do not carry a A2075G 23S rRNA gene mutation always corresponded to a susceptible phenotype). In contrast, the PPV (i.e., the accuracy of detecting a mutation and the isolate being resistant to either TEL or CLI) was only 92.9 and 78.6% for CLI and TEL, respectively.

No genes currently known to confer resistance to gentamycin or florfenicol were observed in any of the isolates, resulting in a 100% concordance between genotype and phenotype.

Overall, the screening software used to identify ARGs (Table 1) had the most impact on both the detection of a gene and the type of gene identified, whereas WGS input (i.e., raw reads or draft assemblies) did not appear to affect the genotype. For example, sipprverse, a software which uses raw reads as input, performed similarly to starAMR and GeneSeekr, software which use draft assemblies as input. While all three pipelines detected the same number of tetracycline resistance genes, there were discrepancies in the detection of two beta-lactam resistance genes; two blaOXA genes were not detected by GeneSeekr using SKESA assemblies but were detected in the SPAdes assemblies and by Sipprverse. An investigation into the SKESA assemblies of CJ-MBS1118R and CJ-MBS1595R found that the blaOXA gene had not assembled (CJ-MBS1118R) or had only partially assembled (CJ-MBS1595R).

In contrast, KMA, which uses raw reads as WGS input, detected the fewest isolates carrying tet(O) and blaOXA, with 14 fewer tet(O) genes detected than the other pipelines. The pipelines also differed in the detection of genes associated with resistance to various aminoglycosides, as KMA and SRST2 were the only pipelines to detect the sat4 gene while GeneSeekr detected four additional isolates carrying the aadE gene (Table 3). The failure to detect a sat4 gene by sipprverse was due at least in part to the ARG not found in the ResFinder database. The association between phenotypic resistance and the detection of genes for resistance to beta-lactam or aminoglycosides other than gentamycin was not verified in this study. While all of the pipelines detected tetracycline and beta-lactam resistance [i.e., tet(O) and blaOXA, respectively], the ability to distinguish between alleles also varied among the pipelines (Table 3). For example, not all pipelines reported a predicted class of blaOXA, whereas KMA and SRST2 were the only pipelines that did not identify the mosaic allele tet(O/32/O), which was found in four isolates. With the exception of two blaOXA genes, assemblies generated using SKESA and SPAdes reported the same genotype when screened with the same software (i.e., GeneSeekr).

The impact of genome coverage on the ability to reliably detect ARGs from raw reads was evaluated using sipprverse, KMA and SRST2. Detection of an ARG was found to be dependent on the screening software, the ARG and the level of genome coverage (Figure 1). Although KMA and SRST2 could detect blaOXA genes at lower coverage levels compared to sipprverse, all three pipelines failed to detect the gene 100% of the time for a proportion of isolates using 30X average genome coverage. For example, SRST2 was able to detect the gene in >90% of the subsamples for 89.5% of the isolates (n = 194/217) at 30X genome coverage and the gene was missed in only one isolate (0.5%) where it could not be detected in at least 50% of the subsample. In contrast, at 30X genome coverage, the blaOXA gene was detected in <50% of the replicates in 3.8% (n = 8/217) and 5.5% (n = 12/217) of isolates using KMA and sipprverse, respectively (Figure 2).

In comparison to blaOXA, a tet(O) or tet(O/32/O) gene was not only more likely to be detected but it could be detected using significantly lower genome coverage (Figure 1). Of the three software, sipprverse correctly detected a tet(O) gene determinant in the most isolates, using the least amount of genome coverage. For example, there was a 95% likelihood of detecting a tet(O) gene at only 5X coverage, compared to a likelihood of 60 and 63% using KMA and SRST2, respectively (Figure 1). For example, using KMA false positives were identified in an increasing number of subsamples as genome coverages increased for three isolates, i.e., gene detected in the subsample but not using all available reads. Although the presence of a tet(O) gene was not detected using KMA and all available reads, all three isolates were phenotypically tetracycline resistant and a tet(O) gene was detected using other pipelines. Using sipprverse, one isolate was reported as a false positive in a proportion of the subsampled reads, identified in up to 20% of the subsamples between one times and 20X coverage. The false positive was not detected by any other pipeline (Table 1) and further investigation using both the assemblies and raw reads suggests the gene is non-functional. The result of a BLAST search of the gene indicates this was likely due to the insertion of a transposon, which corresponded with the isolate being phenotypically susceptible to TET.