Four C. hepaticus strains were isolated from organic pasture raised laying hens in Georgia, United States, in 2021 (Becerra et al., 2023b). Bile samples were collected immediately post-mortem from different birds associated with the outbreak farm and showed the presence of spotty livers on necropsy. Volumes of the bile (30 μL) were plated directly on blood agar plates (TSA + 5% Sheep blood, Remel) and cultured under microaerophilic conditions (Mitsubishi, MicroAero Gas generator, Remel) in an anaerobic jar at 37°C for 7 days. Suspect colonies were picked and suspended in sterile water to make DNA using the boil prep method (10 min at 100°C, followed by cooling and centrifugation to extract the DNA). Confirmation was carried out by PCR targeting the glycerol kinase gene using primers and conditions previously described by Van and collaborators (Van T. T. H. et al., 2017).

Genomic DNA was extracted from each strain using the MagAttract HMW DNA Kit (Qiagen, Hilden, Germany). The DNA sample was quantified using Qubit^® dsDNA HS Assay Kit (Life Technologies, Carlsbad, CA), and the quality was evaluated using a NANODROP 2000 (ThermoFisher Scientific, Waltham, MA), followed by analysis using E-Gel SizeSelect 2% Agarose Gel (Invitrogen, Waltham, MA). The SMRTbell Express Template Prep Kit 2.0 (Pacific Biosciences, Menlo Park, CA) was used to prepare the genomic library for PacBio sequencing with BluePippin (Sage Science, Beverly, MA) size selection to target a library size of >7 kb. Genomic sequencing was performed on the PacBio Sequel II (Pacific Biosciences).

Sequencing yielded 317,377 PacBio reads totaling 300,011,149 bp (200x coverage), and their quality was examined using FastQC v.0.11.9. ¹ The PacBio reads were assembled with Canu v.2.1.1 (Koren et al., 2017) using the “pacbio-hifi” option, and the assembly produced five contigs with an N50 value of 1,528,390 bases and a GC content of 28.04%. Circlator (Hunt et al., 2015) was used to trim and circularize the contigs, which were further polished with Pilon, v.1.20 (Walker et al., 2014), eventually resulting in one chromosome when possible. The genome sequences were annotated using NCBI Prokaryotic Annotation Pipeline and RAST (Aziz et al., 2008). PlasmidFinder (Carattoli and Hasman, 2020) was used to identify plasmids in the genomes sequences.

Genomes from other Campylobacter species were used to perform genome comparison with the four C. hepaticus genomes isolated in this study. 32 genomes of Campylobacter spp. were used including (seven genomes of Campylobacter coli, three C. fetus, eight Campylobacter jejuni, one C. bilis, and 13 other C. hepaticus). Information regarding strains, host, location, and accession numbers are detailed in Table 1.

The virulence genes were identified by searching directly against Virulence Factor Database 2022 using VFanalyzer (Liu et al., 2022), as well as using the ABRicate v.1.0.1, ² and SEED v.2.0 (Overbeek et al., 2005), available within RAST pipeline. On the ABRicate analysis, the presence of a gene was considered only with >90 and > 80% of coverage and identity, respectively. Antimicrobial resistance (AMR) genes and metal-resistance genes (MRG) were also identified through SEED, which classified the ARG and MRG in the sub-system “Resistance to antibiotic and toxic compounds.” In addition, MEGARes v.3.0 (Bonin et al., 2023) and ResFinder v.4.0 (Bortolaia et al., 2020) were also used to identify the AMR genes. Genome visualization and comparison were performed with GView. ³

For the pangenome analysis, all four C. hepaticus genomes were annotated with Prokka v.1.14.6 (Seemann, 2014), and the output files were submitted to Roary v.3.13.0 (Page et al., 2015), which provides information about the core and accessory genes. A phylogenetic analysis based on the concatenated alignment of SNPs from all Campylobacter genomes mentioned above (Table 1) was performed with CSI Phylogeny v.1.4 (Kaas et al., 2014). CSI Phylogeny uses Burrows-Wheeles Aligner (BWA), SAMtools, and BEDtools to align nucleotide sequences followed by SNP calling, MUMmer to assess genome similarities, and FastTree 2 for maximum likelihood tree estimation (Ahrenfeldt et al., 2017). All default parameters of the software were applied. The final phylogenetic tree was created with the Interactive Tree of Life (iTOL) v.6.5.8 (Letunic and Bork, 2021).

An additional pangenome analysis for Campylobacter was performed to evaluate the association of genes with each Campylobacter species. All 34 Campylobacter genomes mentioned above (Table 1) were annotated with Prokka v.1.14.6 (Seemann, 2014), and the output files were submitted to Roary v.3.13.0 (Page et al., 2015). The gene presence or absence table generated by Roary was used as input file in Scoary v.1.6.16 (Brynildsrud et al., 2016), with a p-value cut-off of 0.05, adjusted via the Bonferroni correction, to establish gene clusters associated with each trait (Campylobacter sp.).

Statistical analysis regarding the genetic differences between C. hepaticus and other Campylobacter species were performed using GraphPad Prism 9. The Kruskal-Wallis test was performed to compare each subsystem parameter between the Campylobacter species.

One C. hepaticus genome was closed into a chromosome (RBCL71delta), and the other three genomes were assembled into contigs (Becerra et al., 2023a). The genome from C. hepaticus RBCL71delta has a GC content of 28.03% and consists of 1,516,079 base pairs. It contains 1,488 predicted protein-coding sequences (CDS) encoding 54 RNAs (9 rRNA genes and 42 tRNAs and 3ncRNAs). The other three C. hepaticus genomes have an average GC content and size of 28.05% and 1,558,574 base pairs, respectively (Table 2). No plasmids were identified in any of the sequenced C. hepaticus genomes.

Using the virulence factor database (VFDB) from ABRicate, 70 virulence genes were identified in the four sequenced genomes. The VFanalyzer from Virulence Factor Database website had predicted 101 virulence genes in the four sequenced genomes, of which an average of 35 were unidentified, classified by the same tool only by their predicted function (Figure 1). Most of those genes are related to adherence, colonization, immune evasion, motility, and glycosylation. The genes identified through both methods are included in Supplementary Table S1. Only 38 genes were identified by both ABRicate and VFanalyzer. Use of two search tools gave better results than a single search tool alone.

The Campylobacter invasion antigen B (ciaB) was identified in all four C. hepaticus isolates, confirming that the isolates belong to this Campylobacter species, since this gene is an essential virulence factor for Campylobacter that plays a key role in the invasion of C. jejuni, C. coli, Campylobacter lari, and Campylobacter rectus into human intestinal cells (Konkel et al., 1999; Fouts et al., 2005; Onozato et al., 2009; Lagier and Threadgill, 2014) and ciaB functions as a target of the type III secretion system (Konkel et al., 1999). In the present study, the ciaB gene identified in the four C. hepaticus RBCL strains was submitted to BLAST and showed a match only with the other three C. hepaticus genomes in the NCBI database (C. hepaticus HV10, C. hepaticus UF2019SK1, and C. hepaticus USA52) with an identity and query cover of 100%.

Using the SEED Viewer from RAST Annotations, five genes related to bacteriocin production were identified, including the cvpA gene that encodes the colicin V production protein (Table 3). Genes related to multidrug efflux pumps were also identified, including the cmeABC operon, a system that has been extensively studied in C. jejuni species (Lin et al., 2002, 2003). The four sequenced genomes also harbored genes that promote resistance to copper, cobalt, zinc, and cadmium, and fluoroquinolones.

A comparison of the genome sequences of the 17 C. hepaticus genomes showed a high conservation profile with an average identity of 99.91% against the reference genome, C. hepaticus HV10 (Supplementary Table S2). Pan-genome analysis identified 2,102 genes classified into the core, shell, and cloud genes (Supplementary Figures S1, S2). The C. hepaticus core was composed of 1,328 genes, while 366 and 408 genes were considered as shell and cloud genes, respectively.

For better comprehension, circular comparison visualization of 9 representative C. hepaticus complete genomes was performed based on a BLAST analysis (Figure 2), followed by a pan genome analysis (Supplementary Figure S3), including the reference genome, C. hepaticus HV10. This analysis found that there were many conserved regions between the 9 representative C. hepaticus complete genomes when compared, including some identical gaps. However, next to the 1,000 kbp region, the four genomes sequenced in this study showed a distinct gap compared to the reference genome. This same gap was also observed in the C. hepaticus S11-0036, from the United Kingdom. The C. hepaticus USA52 also has an exclusive gap next to the 1,300 kbp region, observed only in this genome, likely due to some gene deletion. This gap is associated with 18 genes identified in all other C. hepaticus genomes, most were codified as non-identified proteins, which were classified as hypothetical proteins in the annotated genomes. However, in this region, some genes codified to the LexA family transcriptional regulator, the type II toxin-antitoxin system RelE/ParE family toxin, and the helix-turn-helix domain-containing protein, which are not present in the C. hepaticus USA52 isolate (Supplementary Figure S4).

Using the classification of sub-systems from RAST, a comparative analysis of virulence factors between different Campylobacter species was performed. All 34 Campylobacter genomes were included in this analysis (See Methods section). It was found that compared with other Campylobacter species examined (C. jejuni, C. coli, C. fetus, and C. bilis), C. hepaticus has a significantly lower number of genes classified in the following functional groups: iron acquisition, virulence, disease, and defense (Figure 3). In contrast, the C. hepaticus genomes harbor a significantly greater number of genes related to carbohydrate metabolism than other Campylobacter species, except for C. coli. It was also found that in contrast to all other Campylobacter genomes, the C. bilis VicNov18 genome does not appear to harbor any genes related to the iron acquisition system (Figure 3).

The pangenome analysis using all Campylobacter genomes identified a total of 15.145 genes, classified into soft, shell, and cloud genes. No gene was identified in more than 99% of the isolates to be considered as core genes. The pan-WGAS revealed that a total of 73, 50, and 16 genes were associated with C. hepaticus, C. jejuni, and C. coli, respectively (p < 0.05; Supplementary Table S3). However, the analysis did not reveal any significant associations between C. bilis and C. fetus with specific genes. Although some virulence genes were associated to C. hepticus, such as flgE, flaA, and flab, most of the genes associated were hypothetical proteins, and further analysis linking the Scoary results and this pathogen is required.

Single nucleotide polymorphisms (SNPs) were identified from alignments of the examined Campylobacter whole-genomes. Figure 4 shows that C. hepaticus genomes were grouped into two distinct groups, and C. bilis VicNov18 was located as an intermediate group between them. One cluster consists of 9 C. hepaticus genomes from the UK, while the second cluster includes all other C. hepaticus genomes from the United States and Australia. However, isolate S12-002 from the UK was not grouped in the UK cluster since it had ~1–3 SNPs differences from the UK genomes and no SNPs differences in comparison with the United States and Australia genomes (Supplementary Figure S5). The highest number of SNPs observed among the C. hepaticus genomes was 4, indicating that this species has a well-conserved genome (Supplementary Figure S5). It was also found that the C. fetus genomes (C. fetus subsp. fetus and C. fetus subsp. veneralis) are phylogenetically closer to C. hepaticus (~16–19 SNPs) than C. coli and C. jejuni. High bootstrapping values support these data.

This study found that the four genomes characterized (RBCL 71, 76, 81 and 91 delta) differed slightly from the type strain C. hepaticus HV10. Using the closed genome RBCL71delta for comparison it has a genome size of 1,516,079 bp compared to 1,520,669 bp for HV10, this is also reflected in a slightly lower GC content (28.03%) (compared to 28.2% for HV10) (Van et al., 2019). However, a few more coding sequences were detected in the RBCL71 delta genome (1,488) compared to C. hepaticus HV10 (1,474). The four sequenced genomes are slightly larger than the USA52 (Wu et al., 2021) and similar in size to UF2019SK (Arukha et al., 2021), the only two other US C. hepaticus genomes available.

Of interest, no plasmids were detected in any of the four sequenced genomes analyzed in this study or the other two US genomes (Arukha et al., 2021; Wu et al., 2021), while Van and colleagues (Van et al., 2019) only detected plasmids in five of 14 C. hepaticus isolates from Australia that were analyzed, suggesting that mobile elements such as plasmids may not be a common trait of C. hepaticus.

The four new C. hepaticus genomes sequenced in this study would appear to be relatively conserved with a match level of 99.95% with C. hepaticus HV10. This match is greater than that reported by Van and collaborators (Van et al., 2019) where comparison of HV10 with a collection of 22 C. hepaticus genomes consisting of 13 Australian genomes and 9 UK genomes found an identity match of 95.19%. The four genomes sequenced in this study are a first report of C. hepaticus in US pasture raised poultry.

In this study, phylogenetic analysis was performed based on the concatenated alignment of SNPs, and this analysis provides a SNP matrix from the whole-genome alignment. A low number of SNPs was detected among the C. hepaticus genomes, but level present might be enough to initiate an evolutionary process that distinguishes the UK genomes from the C. hepaticus genomes from the United States and Australia. Moreover, the lowest number of SNPs among the C. hepaticus genomes indicates that this species has a conserved genome aside from geographical differences. In addition, considering that C. hepaticus was described in 2016 (Van et al., 2016) and just a few genomes are currently available, the lower number of SNPs is likely plausible. The four C. hepaticus genomes from Georgia (United States) showed 99.95% sequence identity compared to the reference genome C. hepaticus HV10, which sustains the hypothesis that these genomes are highly conserved.

Virulence gene analysis of the four genomes using two different search databases found significant differences in the number of genes identified as virulence genes. The database ABRicate identified 70 virulence genes among the four sequenced C. hepaticus genomes compared to 101 using the VFanalyzer directly from Virulence Factor Database (VFDB) website. When both search engines were compared, only 38 of the identified virulence genes were common to both search tools. Of interest, even with regular updates, the ABRicate database has the same number of sequences since 2019, when the software was curated with a database from 2016 (Chen et al., 2016), while the newer search database VFanalyzer appears more current (Liu et al., 2022) and was updated in 2022. This suggests that the VFDB online resource is more accurate in finding virulence genes if compared by the number of genes found, however use of both search tools resulted in a better overall outcome in identifying all virulence genes in the four genomes examined. This study highlights the importance of using more than one analysis method to evaluate whole genome sequences, since their different approaches alone could provide false-negative results. To our knowledge, this is the first time that two different search databases were compared for output in terms of virulence genes identified for C. hepaticus.

Among the virulence genes identified, the ciaB gene, associated with Campylobacter invasion was found in all four genomes and showed 100% identity match to C. hepaticus HV10, UF2019SK1, and USA52. Similar results were identified by Onozato et al. (2009) when compared the nucleotide sequence of ciaB gene from C. lari isolates with the sequences from C. coli and C. jejuni isolates. The authors suggested that the nucleotide sequence information of the ciaB gene has molecular discrimination efficacy among the Campylobacter species analyzed (Onozato et al., 2009). This finding raises the hypothesis that CiaB might be a well conserved gene among C. hepaticus species and further studies are required to identify its potential as a target for antimicrobial or diagnostics tools.

Although no plasmid or mobile elements were detected in our C. hepaticus genomes, some genes first characterized in C. jejuni and C. coli were identified in the C. hepaticus RBCL genomes. The multidrug efflux pump CmeABC was identified in all four genomes and has also been described by Van and colleagues (Van et al., 2019) in C. hepaticus HV10. This operon is also present in the C. hepaticus strains USA52 and UF2019SK1, from Iowa and Florida, respectively. CmeABC plays an active role in the intrinsic resistance of C. jejuni (Lin et al., 2002), has been extensively characterized in C. jejuni (Lin et al., 2002, 2003; Lin et al., 2005a,b) and has also been detected in C. coli and other Campylobacter species (Olah et al., 2006; Guo et al., 2010) appearing to have a similar genomic organization in all species. Functional analysis identified CmeABC has a role in the efflux of drugs and other antimicrobial agents out of the bacterial host and contributes to the colonization of the poultry host (Lin et al., 2003). Lin and colleagues (Lin et al., 2005b) demonstrated that in the presence of bile salts, there was elevated expression of cmeABC genes in C. jejuni and increased resistance to multiple antimicrobials. These studies also found bile resistance was a natural function of CmeABC. In the present study, this operon was observed in all Campylobacter species, as well as its repressor gene, cmeR. C. hepaticus and C. bilis have been detected and recovered from bile and liver samples of birds affected with SLD (Van T. T. H. et al., 2017; Phung et al., 2022; Van et al., 2022). The presence of CmeABC in C. hepaticus and C. bilis might be an essential factor in these pathogens surviving and colonizing the liver and gall bladder of birds.

In the current study, detection of a full CmeABC operon in C. hepaticus likely plays a key role in the survival and host adaptation of C. hepaticus to the chicken host liver and gall bladder. All four C. hepaticus strains sequenced were isolated from bile samples (Becerra et al., 2023b) and C. hepaticus has been detected in the liver of infected birds with or without characteristic spots using PCR analysis (Becerra et al., 2023b). Given that the CmeABC operon functions to promote resistance to bile then C. hepaticus would appear to have acquired the ability to survive and flourish in this host niche environment.

In concert with CmeABC, genomic analysis also identified an overall reduction in the size of the C. hepaticus genome compared to other members of the Campylobacteriaceae. Of note, reductions were observed in the number of genes associated with the iron acquisition systems (Figure 3). A similar observation was noted for C. hepaticus HV10 (Van et al., 2019). Petrovska and colleagues (Petrovska et al., 2017a,b) also noted a genome reduction in C. hepaticus genomes when analyzed for their association with the poultry host. This loss of system function may be related to the niche environment occupied by C. hepaticus which is the poultry liver and ease of access to specific factors. The liver itself functions to maintain systemic iron balance for the host and also acts as a storage site for excess iron. Therefore this niche may be the perfect environment for C. hepaticus and working in concert with CmeABC may serve as mechanisms to promote and maintain the C. hepaticus population in the liver/bile of infected birds, an observation also supported by Van and colleagues (Van et al., 2019). Further work to validate this observation is required, however, recent work using an iron depletion assay demonstrated the inability of C. hepaticus to grow in the absence of available free iron supporting the observation that the iron metabolism pathways become eliminated when the metabolites are readily available (Petrovska et al., 2017a,b). Also of note, when the C. bilis genome was included in the current analysis it also failed to harbor any iron acquisition systems supporting its role as another C. hepaticus-like strain with a specific niche in the poultry host.

When all four genomes were examined for other metabolism systems, genes associated with carbohydrate metabolism systems were found to have a significantly greater prevalence in C. hepaticus compared to all other members of the genus (C. jejuni, C. fetus and C. bilis) except C. coli. Carbohydrate systems function as pathways that the organism can use to function in a carbohydrate rich environment such as the liver (Petrovska et al., 2017a,b; Van et al., 2019). Generally, these systems are not as common in the genus Campylobacter as the organism is characterized as non-glycolytic (Van et al., 2019) but the presence of greater numbers of genes linked with carbohydrate metabolism may be necessary to cope in a carbohydrate rich environment such as the liver (Petrovska et al., 2017a,b).

Although none of our four sequenced genomes were found to harbor plasmids or the tetO gene which is often found associated with Campylobacter species, the four sequenced strains did, however, show the presence of genes associated with resistance to heavy metals including copper, cobalt, zinc and cadmium and the presence of genes linked with fluoroquinolone resistance (Table 2). Additionally, no genes related to tetracycline or any other antimicrobial class were found in the four C. hepaticus genomes.

One potential reason for the lack of tetO and other antimicrobial resistance genes among the four sequenced isolates could be related to the farm source: all strains originated on an organic farm where no antibiotics are used and if this farm has been in operation for a number of years there may be a lack of persistence of tetO in the absence of a selective pressure. However, some researchers suggest that resistance can still persist in the absence of selective pressure. In the case of the current farm however, little information is available as to the age of the farm and length of time operational as an organic producer but it has been estimated to be more than 5 years, and any comments on resistance persistence would be speculation.

The present study is the first to perform a whole genome comparison among C. hepaticus genomes. The comparison between C. hepaticus and other Campylobacter sp. genomes showed that the lowest number of SNPs was observed between C. hepaticus and C. fetus genomes. This close phylogenetic relationship was never identified before, probably because previous studies examining C. hepaticus phylogeny failed to include C. fetus in genomes their analysis. Some studies included C. fetus in their phylogenetic analysis (Petrovska et al., 2017b; Gregory et al., 2018); however, their analyses were based on 16S rRNA genes and rMLST, respectively, instead of the whole genome. In the current study, C. fetus was included because of its ability to colonize humans and livestock, and is a potential zoonotic pathogen (Chala et al., 2021). Iraola and colleagues (Iraola et al., 2017) demonstrated that C. fetus may have originated as a human pathobiont around 10,500 years ago and may have “jumped” into cattle and adapted as a venereal and abortive pathogen during the livestock domestication period.

Compared with all other Campylobacter genomes, C. hepaticus, C. bilis, and C. fetus have the lowest number of genes associated with iron acquisition and metabolism, while C. jejuni and C. coli have a similar number of genes for this system. C. fetus infects the intestinal tract of several mammalian species and induces infertility and abortion in cattle and sheep due to its ability to colonize and translocate to the placenta (Grogono-Thomas et al., 2000; Iraola et al., 2016). Chemotaxis assays showed that in high concentrations, C. fetus is a chemoattracted to ferrous iron, pyruvate, fumarate, L-serine, L-aspartate, and L-glutamate, but hormones such as estradiol and progesterone are non-chemotactic for this species (Haas et al., 2021). The placenta is an environment rich in iron since this is essential for development of the fetus (Peter, 2013; Cao and Fleming, 2016), and its composition might support the translocation of C. fetus from the gastrointestinal tract to the placenta by chemotaxis. Since this environment is a free nutrient source, the bacteria do not need a complete apparatus for iron metabolism, which explains the lower number of genes related to iron acquisition and metabolism as observed in this study.

In addition, C. fetus also harbors genes associated to resistance to fluoroquinolones, copper homeostasis and tolerance, and cobalt-zinc-cadmium resistance (Table 2), genes that were also evident in C. hepaticus. All these similarities could help explain the close phylogenetic relationship between C. hepaticus, C. bilis and C. fetus species and their ability to hide out in exclusive niches of their animal host.