Campylobacter hepaticus HV10 was published as a draft genome (Van et al., 2017b). In this study, the complete closed genome of HV10 was obtained by combining short Illumina reads, and long PacBio reads and a series of bioinformatics pipelines as described in Lacey et al. (2018). The closed C. hepaticus HV10 genome was deposited in the NCBI database (accession number: CP031611.1). Fourteen Australian C. hepaticus isolates, each from an independent SLD outbreak event, were sequenced (methods as described in Van et al. (2016) and compared with the publicly available whole genome sequences of nine C. hepaticus isolates from the United Kingdom, five C. jejuni, and five C. coli representative genomes extracted from the NCBI database (Table 1). Genomes were assembled using the A5MiSeq pipeline version 20150522 (Coil et al., 2015) and they were annotated using both Prokka 1.14-dev and RAST version 2.0 (Aziz et al., 2008; Seemann, 2014). All assemblies and read sets were deposited in NCBI (Bioproject PRJNA485661).

Homolog identification and pan genome investigation of all predicted coding sequences was performed using two independent methods. Clustering using USEARCH v10.0 at 70% similarity across 70% protein length within the Bacterial Pan genomes analysis (BPGA v1.3) tool pipeline (Chaudhari et al., 2016), and with Roary v3.12 at a protein identity of 70% (-i70) and no splitting of paralogs (-s) (Page et al., 2015).

For phylogenetic inferences, single nucleotide variants (SNVs) were called by aligning reads to a reference genome, C. hepaticus strain HV10, using Snippy v3.2 (https://github.com/tseemann/snippy). Gubbins v2.3.4 (Croucher et al., 2015) was used for the detection and removal of recombinogenic regions, and PHASTER was used to screen for prophage integrations that would be outside the clonal frame (Zhou et al., 2011). A Maximum-likelihood tree was built in RAxML v8.2.12 (Stamatakis, 2014) using the general-time reversible model (GTRCAT) with 1,000 bootstrap replicates. Clustering of strains was performed using RAMI at a patristic distance threshold of 0.05 divergence (Pommier et al., 2009).

The map of the DNA features of C. hepaticus reference strain HV10 was produced in DNAplotter v16.0.0 (Carver et al., 2009). The core and pan genome plots of the 24 C. hepaticus isolates and COG distribution plot of functional categories for coding sequences within these 24 genomes were produced using BPGA.

The C. hepaticus HV10 genome was examined for potential virulence genes by searching against the Virulence Factor Database (http://www.mgc.ac.cn/VFs/main.htm) in ABRicate (https://github.com/tseemann/abricate) (data assessed June 2018) and by inspecting the annotated genome manually. To elucidate the genetic potential of C. hepaticus to cause SLD a pan genome wide association study (PGWAS) was performed using Scoary v1.6.10 (Brynildsrud et al., 2016), and each gene in the C. hepaticus pan genome for association to SLD was screened. Genes of interest were identified as specific to the C. hepaticus genomes (present in 100% of isolates and absent in all reference C. coli and C. jejuni strains). The functionality of these unique genes was inferred from matches to the Pfam database (Finn et al., 2014), Interproscan (Jones et al., 2014), Swiss-Prot (Bairoch and Apweiler, 1997) and Uniprot (UniProt Consortium, 2015). Product descriptions were assigned with homologs of 70% similarity across 90% of protein length. CRISPRFinder (v2017-05-09) (Grissa et al., 2007) was used to analyse CRISPRs.

The annotated genomes of the 23 C. hepaticus isolates were first manually inspected for any potential acquired genetic materials. Contigs with genes annotated as suspected plasmid elements were Blasted against the NCBI database. Significant matches were determined by matches >90% coverage and identity. ABRicate v0.8.7 was used to screen for antibiotic resistance genes.

C. hepaticus HV10 (Van et al., 2016) were grown on Brucella agar (Becton Dickinson) with 5% horse blood (HBA) and incubated at 37°C in microaerobic conditions using CampyGen gas packs (Oxoid).

SLD in chickens was induced by challenge with C. hepaticus HV10. The animal experimentation was approved by the Wildlife and Small Institutions Animal Ethics Committee of the Victorian Department of Economic Development, Jobs, Transport and Resources (approval number 14.16). Hy-Line layer hens (26-weeks old, sourced from a farm that had not observed any SLD in their flocks for several years) were used in the study. Birds were also tested for C. hepaticus to ensure they were C. hepaticus negative before the trial by using the specific PCR developed by Van et al. (2017b) on the cloacal swab samples. Experimental chickens were challenged as previously described (Van et al., 2017a). Briefly, birds were challenged by direct oral gavage with 1x10⁹ CFU of C. hepaticus HV10 strain in 1 ml of Brucella broth, whereas the control chickens were given 1 ml of Brucella broth. The birds were sacrificed 5 days post-challenge and the livers were examined for lesions. Bile samples from all chickens were taken aseptically from the gall bladder and placed in tubes containing RNAlater (Qiagen) for RNA isolation. Samples were kept on ice, transported to the laboratory and processed immediately.

For RNA isolation from C. hepaticus grown in HBA (in vitro), C. hepaticus was harvested from HBA plates and resuspended in Brucella broth to an OD₆₀₀ of 0.5 then centrifuged. The cell pellet was resuspended in RNAlater to stabilize RNA. RNA was extracted using the ScriptSeq Complete kit (Epicenter) following the manufacturer's instructions. The RNA was treated with DNase I (NEB) to remove DNA contamination. The quality of the total RNA in the samples was checked using a Nanodrop spectrophotometer (Thermofisher). RNA was also electrophoresed on a 1% agarose gel and PCR amplified using SLD specific primers (to check genomic DNA contamination) as previously described (Van et al., 2017b). RNA concentration was measured using the Qubit RNA Assay Kit (Life Technologies). The RNA samples were stored at −80°C. Both in vitro and in vivo samples were done in triplicate.

Ribosomal RNA was first removed from the RNA samples using a Ribo-Zero Magnetic Kit (Bacteria) (Illumina). Libraries for Illumina sequencing were prepared using a ScriptSeq v2 RNA-Seq Library Preparation Kit (Illumina) from the rRNA-deleted RNA. The libraries were sequenced on an Illumina MiSeq platform using 300 bp paired end reads.

The Illumina reads were mapped to the reference genome and differentially expressed genes (DEGs) were identified. Raw reads were quality trimmed using Trimmomatic version 0.36 (Bolger et al., 2014), and the trimmed reads were aligned against the C. hepaticus HV10 reference genome using BWA (Li and Durbin, 2010). The SAM files were imported into Blast2Go version 4.1.9 for unique read counts and differential expression analysis (Conesa and Götz, 2008). Parameters for classifying significantly expressed genes (DEGs) were ≥2-fold differences in the transcript abundance and ≤ 0.5% false discovery rate (FDR). The list of up-regulated/down-regulated genes were motif scanned to investigate their biological significances and SEED viewer was used for subsystem functional categorization of the predicted open reading frames (ORFs) from RAST annotation. DEGs were further examined by determining the KEGG Biosynthesis pathway to which they belonged.

Campylobacter jejuni strain 81116 (NCTC11828), C. coli NCTC 11366 and three C. hepaticus isolates, C. hepticus HV10, C. hepaticus 19L and C. hepaticus 44L were used in glucose utilization studies. After cultures were grown in HBA for 3 days, cells were collected and resuspended in physiological saline (0.9% NaCl) to an OD₆₀₀ of 1.0. The medium used to test the ability of C. hepaticus to utilize glucose consisted of inorganic salts (IS) as described previously (Alazzam et al., 2011). L-cysteine (0.2 mM) was used as a nitrogen source, and α-D-glucose (10 mM) (Sigma) was used as the sole carbon source. The experiment was carried out in 24-well plates. Each well-contained 100 μl of culture (OD₆₀₀ = 1). Controls included culture in IS plus L-cysteine only and IS plus α-D-glucose only. 2,3,5 tetrazolium chloride (TTC) (Sigma) (0.0665 g/L) was used as an indicator in all wells (Menolasino, 1959). A color change is observed in growing cultures, indicating utilization. The color change results were read after 36 h of incubation at 37°C. C. hepaticus, C. jejuni, and C. coli were also grown in Brucella broth to confirm their viability. The experiment was repeated twice, each time in biological triplicate.

By combining short Illumina reads and long PacBio reads, the complete closed genome of C. hepaticus HV10 strain was obtained in this study (NCBI accession number CP031611). The genome size of C. hepaticus HV10 is 1,520,669 bp with a GC content of 28.2%. It has been used as the reference genome to compare with the other draft genome sequences.

Genome sizes of the C. hepaticus isolates ranged from 1.475 to 1.565 Mb whereas the representative C. jejuni genomes ranged between 1.617 and 1.800 Mb and C. coli genome sizes ranged between 1.668 and 2.034 Mb. There were between 1472 and 1595 annotated protein-coding genes predicted to be encoded by the 23 C. hepaticus isolates, whereas that were 1622–1954 for five C. jejuni isolates and 1672–2155 for five C. coli isolates (Table 1). A total of 1,059 core genes were conserved in all 33 genomes (fourteen Australian C. hepaticus isolates, nine United Kingdom C. hepaticus isolates, five C. jejuni and five C. coli). Maximum likelihood (ML) phylogeny as produced from RAxML was inferred from 33,157 SNVs. The core genome tree of the 33 genomes comprised 3 phylogenetically distinct lineages corresponding to each of the 3 species (Figure 1A).

Within each phylogroup the mean pairwise nucleotide divergence between genomes was ~3.58%, but the C. hepaticus genomes showed much less divergence (0.46%) than the C. coli (6%) and C. jejuni (4.29%) genomes. Nucleotide divergence from C. hepaticus to C. coli and C. jejuni was ~68 and 63.8%, respectively (calculated across core sequence alignments). The nucleotide divergence between C. coli and C. jejuni was 23.16%, indicating that C. coli and C. jejuni are more closely related to each other than to C. hepaticus. Gubbins did not detect any recombinogenic regions, showing there is no evidence of recombination between the phylogroups, adding further support to the clear separation of C. hepaticus from C. jejuni and C. coli. Accessory genome variation also supports the separation of C. hepaticus from the other species, demonstrating that each phylogroup is a discrete bacterial population that is evolving independently, with limited homologous recombination between groups.

Comparison of the genome sequences of the 23 C. hepaticus isolates revealed a total of 1,360 core genes conserved across all genomes. Maximum Likelihood (ML) phylogeny was produced from 4,812 core genome SNVs and revealed a shallow branching population structure with high bootstrapping support. There is a high level of conservation within the genomes (median of 95.19% coverage of the reference strain HV10), possibly due to the specific niche adaption of the species. The C. hepaticus genomes clustered into five phylogenetic lineages based on the core genome ML tree using a patristic distance of 0.05 (Figure 1B).

The Australian isolates formed two lineages, which differed by on average ~4,429 SNPs. UK isolates formed three lineages. The Australian HV10 phylogroup differ from the UK isolates by between ~500 and 1300 SNPs. The UK lineages, represented by UK cluster 1 and UK cluster 2, differ by ~1,100 SNPs while UK lineage 3 (single isolate) differs to the rest of the UK isolates by ~1,100 to ~1,300 SNPs. Within each phylogroup the SNP frequency was very low with an average of ~41 SNPs (range 0–115 SNPs).

A pan genome of 1,709 unique protein-coding sequences was identified across the 23 C. hepaticus isolates. A map of C. hepaticus HV10 DNA features is presented in Figure 2A. The rapid plateauing of the gene accumulation curve (Figure 2B) revealed an almost closed pan genome, suggesting most of the genetic diversity has been discovered, despite the relatively small sample size of sequenced genomes. Each genome carried on average 103 accessory genes (range = 70–192), and most of those genes were associated with multiple strains, with very few rare/unique (single isolate) accessory genes. This supports the hypothesis of genome reduction/speciation events in C. hepaticus as suggested by Petrovska et al. (2017). Isolates with >100 accessory genes were found to have sequences associated with mobilizable tetracycline resistance plasmids, highly similar to Campylobacter plasmids pCJDM210L and pCC31. The COG distribution plot of functional categories for coding sequences within the 24 C. hepaticus genomes (Figure 2C) showed that transport and metabolism categories of various substrates are mostly encoded in the core genome, while most of the accessory and unique genome variation is categorized as replication, signal transduction, cell wall biogenesis, and motility.

Initial annotation and characterization of the C. hepaticus isolates, using the type strain HV10 as the reference sequence for the species, showed the typical structure for a Campylobacter genome with many genes encoding chemotaxis (11 genes), motility (47 genes), adherence/surface protein (59 genes), as well as various metabolism loci for acquisition of metals and carbohydrates (Table 2).

Various methods of genome comparison have been used to investigate coding sequences that are unique to C. hepaticus. The strongest associations with high specificity and selectivity to C. hepaticus were genes with predicted roles in chemotaxis, capsule and lipooligosaccharide synthesis and metabolism. Four chemotaxis proteins with < 88% homology to known chemotaxis proteins were characterized, which could play a role in the movement of C. hepaticus from the gastrointestinal tract to the liver (Table 2). Significant variation was also characterized in the lipooligosaccharide locus (LOS), a region of the Campylobacter chromosome known to undergo rearrangements and recombination events (Parker et al., 2005; Revez and Hänninen, 2012). Two points of interest in this locus were exclusive to C. hepaticus. Firstly the ganglioside mimics (NeuABC) are rearranged outside of the locus as normally seen in C. jejuni (no longer located between the Waac/WaaM to WaaV/WaaF). Secondly, there was an apparent ~6.6 kb insertion of seven CDS into the cst-II gene, all with functions predicted as various glycosyltransferases. This insertion in the middle of the locus resulted in the truncation of cgtA (Table 2). Roughly 2 kb of the inserted sequence is unique to C. hepaticus, with the remaining 4.6 kb showing high sequence divergence to C. jejuni isolates.

A glucose utilization operon was found to be associated with SLD and is discussed in detail in a later section. All C. hepaticus isolates encode a region of CRISPR-cas genes (type II-cas9 CRISPR), however a CRISPR array (section of repeats and spacers) was only found (CRISPR-finder) in two isolates from the divergent Australian Cluster 2 (three direct repeats and 2 spacers, isolates 19L and 54L). The remaining 22 isolates did not encode a complete CRISPR array, just the cas genes (cas9, cas1 and a fragmented cas2). A region of 13 kb is inserted within the two cas2 CDS of the 22 remaining isolates; ~7kb is unique to C. hepaticus. The GC content of this region is similar to that of the rest of the chromosome (28.03%), with most genes (11 CDS) having unknown functions. Within this region, three CDS encoding for luxA repressor, XRE family transcriptional regulation and type II toxin-antitoxin system mRNA interferase are present. A screen for prophage using PHASTER did not identify any complete prophage integrations within any of the genomes.

Plasmids are present in five out of fourteen C. hepaticus Australian isolates (Table 3). Using ABRicate to screen for antibiotic resistance genes, a single antibiotic resistance gene, tetO, was found in eight isolates (5 from Australia and 3 from the UK), which correlated directly to the presence of plasmid elements. Distinct plasmids were found based on the country of origin of the isolates. UK isolates contained a plasmid highly homologous to the previously characterized C. coli plasmid pCC31 (99% coverage and identity), while the Australian isolates contain plasmids homologous to the C. jejuni pCJDM210L plasmid (93% coverage and 99% identity). This plasmid harbored a type IV secretion system along with a tetracycline-resistant gene. Five of the Australian isolates within this study (27L, 84B, Ace1, Ace8659, and AceM3a) carry the plasmid and it accounts for roughly half of the gene content within the accessory genome of these isolates. As short-read sequence data was used it was not possible to assemble the plasmid in its entirety. At least three contigs from each of these genomes were highly conserved and carried plasmid elements.

C. hepaticus encodes three ribosomal RNA operons, however two have been disrupted by the insertion of multiple CDS between the 16S and the 23S genes. A glucose utilization operon and an oligopeptide transporter operon were located within ribosomal RNA operons (Figure 3). Analysis of the insertions showed that the glucose utilization and oligopeptide transporter regions have GC content of 28.17 and 27.56% respectively, which is similar to the average GC content of the HV10 genome, 28.2%.

C. hepaticus cultures incubated for 24 h in IS media containing L-cysteine and D-glucose-6-phosphate, or L-cysteine and D-glucose showed color development (due to TTC) and therefore demonstrated utilization of the substrates, whereas C. jejuni and C. coli cultures did not. There was no color development in C. hepaticus cultures incubated in IS plus L-cysteine only or IS plus D-glucose or D-glucose-6–phosphate. This was due to a lack of carbon source and nitrogen source, respectively. The color change was observed in all cultures grown in Brucella broth, demonstrating the viability of the inoculated cultures.

In C. hepaticus recovered from the gall bladder of SLD experimentally infected birds, 410 genes were differentially expressed (False Discovery Rate (FDR) < 0.05) when compared to in vitro grown bacteria. There were 164 up-regulated genes ((log2-fold-changes) > 1.0) in vivo and 246 down-regulated genes (logFC < −1.0). Functional gene categorization assessed using the SEED Viewer, showed that the 410 differentially expressed genes belonged to 56 subcategories (Figure 4). Notably, all genes associated with polyhydroxybutyrate (PHB) metabolism (Figure 5) were up-regulated (EC 1.1.1.30: D-beta-hydroxybutyrate, EC 2.3.1.9: Acetyl-CoA acetyltransferase, EC 2.8.3.5: Succinyl-CoA:3-ketoacid-coenzyme A transferase, and genes encoding D-beta-hydroxybutyrate permease, short chain fatty acids transporter and 3-ketoacyl-CoA thiolase). These genes may play a role in stress response in C. hepaticus and are putative virulence factors (Table 2).

The gene clusters encoding Ni-Fe-hydrogenase were up-regulated in the cells recovered from bile (Table 2). Six of out eight genes associated with nitrate and nitrite ammonification (nitrogen metabolism system) were also up-regulated and only one was down-regulated in bile samples (Figure 4B). Transcripts from the phosphate transport system of pstS and pstC were increased in abundance in vivo compared to in vitro (Table 2). RNA-Seq identified increased abundance of many transcripts associated with copper homeostasis and up-regulation of pathogenesis-associated glutamine ABC transporters, papP and papQ (Table 2). The neuB (N-acetylneuraminate synthase) and neuC (UDP-N-acetylglucosamine) genes, necessary for sialic acid synthesis, were both up-regulated in the bile-derived bacteria. Thirteen genes associated with flagella motility were down-regulated, and only five were up-regulated. Increased expression of flagella associated genes included genes encoding flagella motor protein (MotA) and flagella switch motor protein (FliN). Down-regulated genes included a putative lipoprotein required for motility, a motility integral membrane protein and flagella-associated genes including flgB, flgD, flgF, flgG, flgH, and flgI. In addition, many genes in the aromatic amino acids and derivatives category were down-regulated (Figure 4B), including genes involved in common pathways for synthesis of aromatic compounds, tryptophan synthesis, and chorismate synthesis (intermediate for synthesis of tryptophan) and none of the genes in this category were up-regulated. Similarly, genes associated with the production of methionine (lysine, threonine, methionine, and cysteine subcategory, Figure 4B) were down-regulated. RNA-Seq identified decreased abundance of transcripts associated with tRNA processing, RNA methylation, and RNA pseudouridine syntheses. On the other hand, there was variation in the expression of genes involved in the oxidative phosphorylation pathway, as genes encoding the enzyme NADH dehydrogenase (EC 1.6.99.3) were down-regulated, while many genes encoding enzymes NADH ubiquinone oxidoreductase (EC 1.6.5.3) and ubiquinol-cytochrome C reductase (EC 1.10.2.2) were up-regulated (Table 2).

AA and PS were employed by company Scolexia Pty Ltd. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.