BACKGROUND
the pathogen campylobacter fetus  can cause disease in both animals and humans.
in humans, c. fetus infections vary from acute diarrhea to systemic illness  <cit> .
in animals, c. fetus infections can cause abortion and infertility, mainly in cattle and in sheep  <cit> .
the mammal-associated c. fetus subspecies are c. fetus subsp.
fetus  and c. fetus subsp.
venerealis   <cit> , whereas c. fetus subsp.
testudinum  is associated with reptiles  <cit> .
c. fetus subsp.
venerealis includes a biochemical variant, designated c. fetus subsp.
venerealis biovar intermedius   <cit> .
c. fetus subsp.
venerealis has been described to be the causative agent of bovine genital campylobacteriosis , associated with infertility and abortion in cattle  <cit> .
bgc is notifiable to the world organisation for animal health .
a crucial element in the bgc control program relies on the subspecies identification of c. fetus isolates.
currently, the methods prescribed by the oie to differentiate cff, cfv and cfvi are tolerance to 1 % glycine and h2s production  <cit> : cff is tolerant to 1 % glycine and able to produce h2s, cfv is not tolerant to 1 % glycine and not able to produce h2s and cfvi is not tolerant to 1 % glycine  and able to produce h2s .
the biochemical tests are hampered by poor reproducibility  <cit> , and the phenotypes are not completely consistent with the genomic characteristics of the c. fetus strains, since phenotypically-identified cff strains were genotypically identical with cfv strains  <cit> .
an obvious distinguishing and important feature of c. fetus cells is the surface layer , which is considered to be associated with the pathogenicity of c. fetus strains  <cit> .
c. fetus cells can express two types of surface array proteins , which correlates with the serotypes of the bacterium; cfv strains are serotype a; cff strains may be either serotype a or serotype b and rarely serotype ab; and cft strains are serotype a, serotype b or serotype ab  <cit> .
the molecular method multi-locus sequence typing  was recommended to differentiate cff and cfv strains  <cit> .
however, a recent study showed that the current mlst scheme was not able to reliably differentiate the c. fetus subspecies, as a cff strain was isolated with the cfv-associated mlst st- <dig> genotype  <cit> .
whole-genome analysis provides fine-scale resolution of bacterial genomes and allows the calculation of evolutionary events, as shown for c. coli and c. jejuni  <cit> .
whole-genome analysis has improved, and will continue to improve our understanding of the features that distinguish c. fetus subspecies and the evolutionary forces that have acted on c. fetus over time.
in this study, we performed a core genome single-nucleotide polymorphisms  analysis of  <dig> c. fetus cff and cfv genomes to identify subspecies-specific snps.
we performed a snp-based phylogenetic analysis of the core genomes and a beast analysis to estimate the divergence dates of cff and cfv strains.
additionally, we investigated whether the genomes contain specific snps or genes that could be associated with the biochemical tests and different clinical features of the c. fetus strains.
methods
bacterial strains and whole genome sequencing
in this study,  <dig> c. fetus strains from different countries and sources were included .
the strains were biochemically characterized , using h2s production in medium amended with  <dig>  % cysteine-hcl and 1 % glycine tolerance, as described before  <cit> .
genotypic subspecies characterisation was performed using mlst  <cit>  and aflp  <cit> .
the whole genome sequence data of all c. fetus genomes was published previously  <cit> .
briefly, the c. fetus genomes  were sequenced using a roche  <dig> gs-flx+ genome sequencer with titanium chemistry and assembled into contigs using the newbler assembler .
the genomes with accession numbers starting with err were sequenced using an illumina hiseq and assembled into contigs with the velvet assembler in an established pipeline at the sanger institute.
the genomes of four c. fetus strains  were closed by using a pacbio rs sequencer and assembled into contigs using quiver  with the base calls validated using illumina miseq reads.table  <dig> strain characteristics
n.a., not available; + positive; − negative; cff, c. fetus subsp.
fetus; cfv, c. fetus subsp.
venerealis; cfvi, c. fetus subsp.
venerealis biovar intermedius; cft, c. fetus subsp.
testudinum
the sequences with accession numbers starting with err are available from the european nucleotide archive   <cit> , and the remaining sequences are available from ncbi genbank .
the following genomes were also included; cff strain h1-uy , cfv strain 642– <dig> , cfv strain b <dig> , cfv strain nctc  <dig>  and cfvi strain 99/ <dig> .
the quality of the whole genome sequence data was assessed with the checkm tool  <cit> , showing a suitable completeness score of >96 % with  <dig> tested marker genes.
genome alignment and phylogenetic core genome snp analysis
for phylogenetic core genome snp analysis, whole genome sequences of  <dig> c. fetus isolates were aligned using parsnp v <dig>   <cit> .
we included cft strain 03- <dig> as an outgroup for the phylogenetic core genome snp analysis, but excluded this strain from visualization of the phylogenetic core genome snp tree to get a better resolution of the cff and cfv branches.
snp discovery was focused on a comparison between cff and cfv strains; therefore, cft strain 03- <dig> was replaced by cff strain 82- <dig> as a reference for this analysis.
recombination regions in the core genome alignment were detected and visualized using gubbins  <cit> .
a phylogenetic tree was constructed using fasttree <dig>  <cit>  with a generalized time-reversible model and gamma correction on the recombination-filtered snps in the core genome of all isolates, including cft strain 03- <dig>  the resulting tree was rooted on cft strain 03- <dig> using dendroscope  <cit>  prior to visualization using itol  <cit> .
beast analysis
recombination-filtered non-synonymous snps of the mammal-associated c. fetus isolates were extracted from the gubbins results and used for divergence dating in beast  <cit> , using the isolation dates as tip dates in the phylogenetic tree, as outlined in the beast manual with the following modifications:  <dig> ,000× sampling and a general time reversible  model plus gamma correction as the distance model.
combinations of a strict clock, log-normal clock, exponential clock and random local clock as the clock model and a constant population, exponentially-growing population, and a bayesian skyline plot with six groups as demographic models were used.
maximum ess values and lowest 95 % confidence intervals  in the predicted divergence dates of the clades were obtained with a random local clock and a bayesian skyline plot with four groups as the demographic model.
calculation of branch-specific dn/ds ratios
to calculate the dn/ds ratio per branch in the phylogenetic tree, we aligned the genomes of the mammal-associated c. fetus isolates with parsnp v <dig>   <cit>  and cft strain 03- <dig> as outgroup.
synonymous and non-synonymous snps were determined on the basis of their location in the coding regions of the cff strain 82– <dig> reference genome.
recombination regions were detected using gubbins  <cit>  and excluded from the alignment.
ancestral state reconstruction of the node sequences was performed using fastml  <cit>  with a generalized time-reversible model and gamma correction.
the dn/ds ratios were determined per branch between node sequences.
comparison of snps and genes with phenotypic characteristics of the strains
as traditional differentiation of c. fetus subspecies are based on the biochemical 1 % glycine tolerance and h2s production tests, the genomes were screened for genes and snps that were associated with these phenotypic characteristics, and tested for 1 % glycine tolerance and h2s production in cysteine-rich medium as described above.
the protein-encoding gene presence and absence were determined using blast-based all vs all comparisons with prokka-annotated genomes  <cit>  using roary  <cit> , which clustered the proteins using mcl-edge  <cit> .
a fisher’s exact test was used to calculate the two-tail probability value  of respectively the detected snps and genes versus the outcome of the biochemical tests.
calculation of clinical association
a fisher’s exact test was used to calculate the two-tail probability value  of snps and genes which were specifically found in strains that were isolated from bovine abortions.
significantly-associated snps or genes were checked for their presence in cff strain 04/ <dig>  as this is a non-cfv bovine clinical strain.
snps and genes that were not present in this strain were excluded from this analysis, to separate clinical-associated from phylogenetically-associated snps and genes.
availability of data
genome sequences are available from the european nucleotide archive  and from ncbi genbank, with the accession numbers listed in table  <dig>
RESULTS
phylogenetic analysis of the core genome snps
phylogenetic analysis, based on core genome snps using cft as the outgroup, showed that the mammal-associated cff and cfv genomes group into five distinct clades .
the division of clades was consistent with the serotypes of the strains; clades  <dig> and  <dig> consist of serotype b strains and clades  <dig>   <dig> and  <dig> consist of serotype a strains.
the division of clades was also consistent with the classification of mlst sequence types , except for cff strain h1-uy, which has the cfv-associated mlst st- <dig> genotype  <cit> .
interestingly, snp phylogeny showed the divergence of clade  <dig> and clade  <dig> from a common cff ancestor.
clade  <dig> consists of phenotypically-identified cff, cfvi and cfv strains  without a clear separation in the phylogeny.
however, when the genotypic characterization is used  <cit> , clade  <dig> consists exclusively of cfv and cfvi strains and clade 1– <dig> of cff strains.fig.
<dig> ml tree of c. fetus core genome and mammalian associated c. fetus recombination regions.
cft 03- <dig> was excluded from the visualization.
tree branches are labeled with the clade numbers, serotype and mlst sts of the c. fetus strains.
subspecies identification , h2s production test results , the presence of the putative cysteine transporter  and the gubbins predicted recombination regions in the mammalian associated c. fetus core genome are shown, with on the right side a zoomed ml tree of clade  <dig>  gene flge <dig> encodes a flagellar hook protein and cjei a type ii r-m system protein.
scale bar at the bottom left represents the nucleotide substitutions per site
beast analysis determined that diversification of the mammal-associated c. fetus is a relatively recent event under disruptive selection pressure
we used beast to determine the divergence date of the different clades .
the topology of the beast tree was slightly different than the maximum likelihood  snps tree , especially for the serotype b strains.
this may be either the result of using only synonymous snps in the beast analysis or because of the differences between ml methods and methods based on coalescence.
the 95 % confidence interval  of the beast analysis was too large to estimate the divergence date of clade  <dig> and clades  <dig>   <dig> and  <dig>  and the divergence date of clades  <dig> and  <dig>  .
interestingly, the beast analysis showed that all currently circulating strains in clades  <dig>   <dig>   <dig> and  <dig> diverged very recently from each other.fig.
<dig> beast tree with divergence dates.
beast analysis 95 % ci divergence dates are shown within brackets and with blue bars on branch nodes.
dn/ds ratios are shown in italics on the branches.
branches are colored according the clock rates, represented as substitutions per site per year.
scale bar at the bottom represents thousand years ago
cfv is restricted to the bovine genital tract, but cff can be isolated from different animal species, e.g.
cattle, sheep, goats, pigs, horses, and humans  <cit> .
therefore, we attempted to investigate whether this change in niche preference has resulted in increased evolutionary pressure on cfv.
we estimated the rate of molecular evolution by comparing the ratio of non-synonymous  versus synonymous  snps .
the dn/ds ratio was the highest in the strains of clade  <dig>  compared to the ratio of the strains of clade  <dig>  and in all branches preceding, based on node sequences reconstructed using fastml .
an increased dn/ds ratio  is indicative for recently diversifying genomes, whereas lower dn/ds ratios  correspond to older snps and more diverged genomes  <cit> .
all clades have a dn/ds ratio >  <dig> , with the highest dn/ds ratio of clade  <dig> of  <dig> , showing that the strains of clade  <dig> are evolving under relatively more diversifying selection than the strains of the other clades.
the estimated molecular clock rate for the mammal-associated c. fetus strains was  <dig> – <dig>  × 10− <dig> substitutions per kb per year, which corresponds with the estimated molecular clock rate of c. jejuni of  <dig> – <dig>  × 10− <dig> per kb/year  <cit> .
clade-specific snps and orthologs
the numbers of snps and orthologous genes  specifically present in each clade  are presented in table  <dig>  as well as the specific snps and orthologs for clades representing the serotypes and subspecies of the strains.
the details and annotation of the detected snps and orthologs are listed in additional file 1: table s <dig> table  <dig> core genome snps specific for c. fetus clades
cff, c. fetus subsp.
fetus; cfv, c. fetus subsp.
venerealis; cfvi, c. fetus subsp.
venerealis biovar intermedius
clades 1– <dig> consisted only of cff serotype b strains.
for clade  <dig>   <dig> snps and  <dig> orthologs were found, specifically present in clade  <dig> strains .
the two clade 1-specific orthologs were hypothetical proteins with unknown function.
specific for clade  <dig>   <dig> snps and  <dig> orthologs were found .
the specific orthologs were located in a region encoding mainly hypothetical proteins, as well as a fic domain protein and a transposase.
this clade 2-specific region is chromosomally located in cff strain 04/ <dig> .
all strains of both clades  <dig> and  <dig>  contain  <dig> specific snps and  <dig> specific orthologs.
the clade  <dig> and 2-specific orthologs included a toxin/antitoxin system , restriction-modification  system-associated genes  and s-layer associated genes .
clades  <dig> and  <dig> consist of cff serotype a strains.
clade  <dig> contains  <dig> specific snps and six specific orthologs .
the specific orthologs encoded hypothetical proteins of unknown function, as well as r-m type i  and r-m type iii genes.
for clade  <dig>   <dig> specific snps and no specific orthologs are found.
all strains of clades 3– <dig>  contain  <dig> snps and  <dig> specific orthologs that are mainly located in one region  that includes a glycosyltransferase, a methyltransferase and a ketoreductase.
clade  <dig> consists of both genotypically-identified cfv and cfvi strains.
for this clade,  <dig> specific snps were found .
the genotypically-identified cfv strains of clade  <dig> contain five specific snps and six specific orthologs.
the specific orthologs encoded hypothetical proteins, a transcriptional regulator , a resolvase and a transposase.
analysis of snp regions that are subject to recombination
interestingly, 42 %  of the snps differing between clades 1– <dig> and clade  <dig> are located in regions <1 kb apart.
this is suggestive of recombination, since on average a snp would be expected to occur approx.
every  <dig>  kb, based on the number of snps identified here in genomes of approx.
<dig>  mbp.
.
we predicted potential recombination regions in the mammal-associated c. fetus core genome alignment with gubbins  <cit>  and visualized the results .
using this analysis, one region showed recombinations in all strains of clades 1– <dig> and three recombination regions were found in all strains of clade  <dig>
the recombination region common to all strains of clades 1– <dig> encoded the flagellar hook gene flge <dig>  .
all strains of clade  <dig> contained a recombination region  with a type ii r-m gene   and genes involved in sugar biosynthesis, including: a nad-dependent epimerase/dehydratase ; a nucleotide sugar dehydrogenase ; and a polysaccharide biosynthesis protein .
another clade 5-specific recombination region is cff8240_1473-cff8240_ <dig> that encodes: leader peptidase a ; ribose-phosphate pyrophosphokinase ; a subunit of aspartate carbamoyltransferase , a formate dehydrogenase subunit  and a glutamate synthase subunit .
in all strains of clade  <dig>  a recombination region was found , encoding: several radical sam domain proteins; two transketolase subunits ; a methyltransferase ; and two glycosyltransferases .
additionally, the s-layer  region  was also identified by gubbins as a recombination region, although the sap locus itself is not included in the core genome alignment due to assembly issues of this region  <cit> .
deletion of a putative cysteine transporter is associated with h2s production-negative strains
currently, the oie-prescribed biochemical tests to differentiate cff, cfv and cfvi consist of determining tolerance to 1 % glycine and production of hydrogen sulfide  from l-cysteine  <cit> .
we analyzed the genomes of the c. fetus strains for the presence of specific genes associated with these phenotypic characteristics.
interestingly, h2s negative strains  did not encode two subunits of an amino-acid abc transporter , that together with an atp-binding subunit  putatively form an amino-acid abc transport system involved in cysteine transport  .
the absence of the two abc transporter genes in the c. fetus h2s negative strains possibly explains the phenotypic characteristic of these strains.
cfvi strain 98/ <dig> was described before as h2s positive  <cit> ; however, in this study, the isolates that were sequenced were biochemically characterized and the sequenced isolate of strain 98/ <dig> was h2s production negative, in accordance with the absence of the putative cysteine transporter in this genome.
conjugative transfer region associated with strains from bovine abortions
the c. fetus strains were arranged according to their clinical features: nine strains were from bovine abortions; seventeen strains were from screenings; and for fifteen strains, the clinical features were unknown .
no snps were found that were both present in the only known cff strain from a bovine abortion  and significantly present in the cfv/cfvi strains from bovine abortions.
however, one region was significantly present in all strains from bovine abortions ; this region contains genes encoding a filamentation induced by cyclic amp  domain protein , a dna-binding protein  and conjugative transfer  proteins .
this region is located in cff strain 04/ <dig>  cfv strain 97/ <dig> and cfv strain 84- <dig> on a megaplasmid or an extra-chromosomal element  <cit> .
this tra/trb region is present in  <dig> c. fetus strains isolated from non-abortion cases; cff strain 98/v <dig> and cfv/cfvi strains ccug  <dig>  642- <dig>  adri  <dig>  b <dig>  84- <dig>  lmg  <dig>  nctc  <dig>  zaf  <dig>  adri  <dig> and 99- <dig> .
discussion
the core genome snp analysis of the mammal-associated c. fetus strains identified a large number of snps.
since available pcr methods and mlst typing schemes are not able to identify c. fetus strains correctly to subspecies level  <cit> , cfv/cfvi-specific snps could be used to develop new diagnostic methods, such as a snp probe-pcr, or to improve the current mlst scheme.
however, two drawbacks of a pcr based on cfv/cfvi-specific snps are that the stability of such snps is unknown and that the design of a reliable pcr based on specific snps requires a high degree of optimization  <cit> .
in a previous study, we presented the division of the mammal-associated c. fetus strains into two different clusters, based on a core gene alignment of five cff and  <dig> cfv/cfvi strains  <cit> .
in this study,  <dig> cff strains and  <dig> cfv/cfvi were included, which gave a better resolution of the branches with cff strains, when compared to a previous analysis with five cff strains  <cit> .
the current core genome snp analysis separated the c. fetus strains into five different clades, that contain either the genotypically-identified cff  or cfv/cfvi  strains.
we demonstrated that the biochemical differentiation of the cff, cfv and cfvi strains is not supported by the core genome snp phylogeny: clade  <dig> consists of phenotypically-identified cff, cfv and cfvi strains without a congruent separation in the phylogeny of these strains this showed that the phenotypic separation of the mammal-associated c. fetus strains is not supported by the core genome snp phylogeny, which gives rise to the consideration if the current phylogenetic subspecies differentiation is still reasonable, as previously mentioned  <cit> .
forces of cff and cfv diversification
phylogenetic analysis indicated that cfv/cfvi clade  <dig> and cff clade  <dig> have a common cff ancestor.
no association of the clades with geographic origin and/or host specificity was observed.
interestingly, the short terminal branches and recent diversification suggest that the isolates within each clade, except for strain h1-uy, have diverged recently.
this is in contrast to what is observed in cft, where terminal branches of the isolates are much longer and a much higher genome diversity is observed  <cit> .
the world-wide spread of the mammal-associated c. fetus strains is potentially associated with the improvements in cattle breeding by selection and cross-breeding which started around the 1700s and 1800s, followed by extensive spread of high-producing dairy cows and beef cattle  <cit> , which may have carried only a very limited number of c. fetus clones.
we suggest that this spread from a limited number of sources has resulted in our current observation of the highly clonal nature of mammal-associated c. fetus.
alternatively, the observation of only recently diverged strains is a result of selective sweeps, when a population member with an advantageous trait will take over the population before the trait can spread to other members  <cit> .
it is unlikely that selective sweeps occurred in the mammal-associated c. fetus population, since cattle are intensively monitored worldwide for the presence of c. fetus subsp.
venerealis and only one non-clonal c. fetus strain was found .
cff strain h1-uy was isolated from the blood of a rural worker in  <dig>  who was diagnosed with cellulitis and was in daily contact with cattle  <cit> .
the rates of changes at non-synonymous and synonymous snps indicate whether a gene is under purifying or diversifying selection  <cit>  and is expressed in the dn/ds ratio.
the dn/ds ratio of bacterial genes under stabilizing selection falls generally within the range  <dig> – <dig>   <cit> .
in the mammal-associated c. fetus, the dn/ds ratios were all >  <dig> , showing that all genomes are under diversifying selection.
the highest dn/ds ratio was found for the clade  <dig> strains  .
this high ratio may be due to the genetic features of these strains: they have no crispr-cas systems and the clade 5-strains have  <dig> specific snps in the type ii r-m system gene cjei compared to the other strains .
type ii r-m systems can undergo major changes in specificity by recombination events as shown for helicobacter pylori  <cit>  and r-m type ii systems may play a role in plasmid transformation, as described for c. jejuni where knockout mutagenesis of gene cjei resulted in a strain with a  <dig> -fold-enhanced transformation efficiency  <cit> .
it is unknown whether these snps influence the functionality or specificity of this type ii r-m system in c. fetus strains, but the clade 5-strains are possibly more susceptible to insertion of foreign dna, since the genomes contain genomic islands  encoding t4sss, phages and insertion sequences  <cit> .
h2s negative phenotype of cfv strains associated with loss of a putative cysteine transporter
the original classification of the mammal-associated c. fetus subspecies cff and cfv are based on the 1 % glycine tolerance test  <cit> .
cfv strains can be discriminated from cfvi strains with the h2s production test  <cit> .
we show that the h2s negative c. fetus strains have lost a putative cysteine transporter.
without the encoded transporter, the cells are possibly less capable of importing cysteine, which under normal situations is reduced while forming h2s.
the h2s-negative cfv strains have a niche restriction to the genital tract of cattle, whereas the h2s-positive cfvi and cff strains are not restricted to the genital tract and assumed to be able to colonize the intestines as well  <cit> .
one may speculate that the cfv strains have a defect causing the incapability to grow outside of the genital tract, but it is unknown if this partial deletion of the putative cysteine transporter can be associated with the niche restriction of these strains.
virulence-associated genes specific for clades and strains from bovine abortions
since the first description of the c. fetus subspecies in  <dig>  it is presumed that cfv cause disease in the genital tract of cows, like enzootic venereal sterility and abortions in pregnant cows, and that cff only cause sporadic abortions  <cit> .
it is unknown which genomic features are responsible for the pathogenicity of the c. fetus strains, but potential candidates are the surface layer of the c. fetus cells  <cit>  and type iv secretion systems   <cit> .
we studied the genomes of nine strains from bovine abortions and observed that, in addition to the s-layer proteins and t4sss, all the strains from bovine abortions contain a region encoding conjugative transfer  proteins.
this tra region is located in the closed genomes of strains 04/ <dig>  97/ <dig> and 84- <dig> on a plasmid/ice, which also contain a trb-t4ss and fic genes.
the fic domain proteins encoded by this region are of potential interest, because fic domain proteins have an immune-modulatory function by influencing the cytoskeleton organization of the host cells  <cit>  and are translocated by t4sss in legionella pneumophila, coxiella burnetii and bartonella henselae  <cit> , which has also been proposed for c. fetus strains  <cit> .
the trb/tra genes are also present in  <dig> c. fetus strains isolated from non-abortion cases .
it is possible that these strains can also cause bovine abortions; however, more studies, specifically animal experiments investigating the virulence of c. fetus strains are required to assess this hypothesis.
CONCLUSIONS
phylogenetic core genome snp analysis divided the mammal-associated c. fetus strains into five different clades, which were consistent with the serotypes, but not with the phenotypes of the strains.
beast analysis showed that the clade with genotypically-identified cfv/cfvi strains has evolved from a cff ancestor under diversifying selection.
phylogenetic analysis of the core genome snps did not differentiate h2s-negative cfv from h2s positive cfvi strains.
the partial deletion of a putative cysteine transporter is observed in all h2s negative cfv strains.
additional files
additional file 1: table s <dig>  table with details of c. fetus snps and genes, specific for clades and bovine abortion strains.
additional file 2: figure s <dig>  schematic representation of the putative cysteine transporter genes cff8240_ <dig> - cff8240_ <dig>  in strain cff  <dig> with flanking genes rard and psga .
the start and stop positions of respectively rard and psga are shown in italics.
cfv strain 97/ <dig> lacks one complete orf and contain a remnant of one orf of the putative cysteine transporter.
shown are the similar deletion sites sequences in the cfv strains containing the incomplete transporter.
abbreviations
bgcbovine genital campylobacteriosis
c. fetuscampylobacter fetus
cffc.
fetus fetus
cftc.
fetus testudinum
cfvc.
fetus venerealis
cfvic.
fetus venerealis biovar intermedius
gigenomic island
kyathousand  years ago
snpsingle nucleotide polymorphism
t4sstype iv secretion system

 the pathogen campylobacter fetus can cause disease in both animals and humans. in humans, c. fetus infections vary from acute diarrhea to systemic illness. in animals, c. fetus infections can cause abortion and infertility, mainly in cattle and in sheep. the mammal-associated c. fetus subspecies are c. fetus subsp. fetus and c. fetus subsp. venerealis, whereas c. fetus subsp. testudinum is associated with reptiles. c. fetus subsp. venerealis includes a biochemical variant, designated c. fetus subsp. venerealis biovar intermedius. c. fetus subsp. venerealis has been described to be the causative agent of bovine genital campylobacteriosis, associated with infertility and abortion in cattle. bgc is notifiable to the world organisation for animal health. a crucial element in the bgc control program relies on the subspecies identification of c. fetus isolates. currently, the methods prescribed by the oie to differentiate cff, cfv and cfvi are tolerance to 1 % glycine and h2s production: cff is tolerant to 1 % glycine and able to produce h2s, cfv is not tolerant to 1 % glycine and not able to produce h2s and cfvi is not tolerant to 1 % glycine and able to produce h2s. the biochemical tests are hampered by poor reproducibility, and the phenotypes are not completely consistent with the genomic characteristics of the c. fetus strains, since phenotypically-identified cff strains were genotypically identical with cfv strains. an obvious distinguishing and important feature of c. fetus cells is the surface layer, which is considered to be associated with the pathogenicity of c. fetus strains. c. fetus cells can express two types of surface array proteins, which correlates with the serotypes of the bacterium; cfv strains are serotype a; cff strains may be either serotype a or serotype b and rarely serotype ab; and cft strains are serotype a, serotype b or serotype ab. the molecular method multi-locus sequence typing was recommended to differentiate cff and cfv strains. however, a recent study showed that the current mlst scheme was not able to reliably differentiate the c. fetus subspecies, as a cff strain was isolated with the cfv-associated mlst st- genotype. whole-genome analysis provides fine-scale resolution of bacterial genomes and allows the calculation of evolutionary events, as shown for c. coli and c. jejuni. whole-genome analysis has improved, and will continue to improve our understanding of the features that distinguish c. fetus subspecies and the evolutionary forces that have acted on c. fetus over time. in this study, we performed a core genome single-nucleotide polymorphisms analysis of c. fetus cff and cfv genomes to identify subspecies-specific snps. we performed a snp-based phylogenetic analysis of the core genomes and a beast analysis to estimate the divergence dates of cff and cfv strains. additionally, we investigated whether the genomes contain specific snps or genes that could be associated with the biochemical tests and different clinical features of the c. fetus strains. methods bacterial strains and whole genome sequencing in this study, c. fetus strains from different countries and sources were included. the strains were biochemically characterized, using h2s production in medium amended with  % cysteine-hcl and 1 % glycine tolerance, as described before. genotypic subspecies characterisation was performed using mlst and aflp. the whole genome sequence data of all c. fetus genomes was published previously. briefly, the c. fetus genomes were sequenced using a roche gs-flx+ genome sequencer with titanium chemistry and assembled into contigs using the newbler assembler. the genomes with accession numbers starting with err were sequenced using an illumina hiseq and assembled into contigs with the velvet assembler in an established pipeline at the sanger institute. the genomes of four c. fetus strains were closed by using a pacbio rs sequencer and assembled into contigs using quiver with the base calls validated using illumina miseq reads.table strain characteristics n.a., not available; + positive; − negative; cff, c. fetus subsp. fetus; cfv, c. fetus subsp. venerealis; cfvi, c. fetus subsp. venerealis biovar intermedius; cft, c. fetus subsp. testudinum the sequences with accession numbers starting with err are available from the european nucleotide archive, and the remaining sequences are available from ncbi genbank. the following genomes were also included; cff strain h1-uy, cfv strain 642–, cfv strain b, cfv strain nctc and cfvi strain 99/. the quality of the whole genome sequence data was assessed with the checkm tool, showing a suitable completeness score of >96 % with tested marker genes. genome alignment and phylogenetic core genome snp analysis for phylogenetic core genome snp analysis, whole genome sequences of c. fetus isolates were aligned using parsnp v. we included cft strain 03- as an outgroup for the phylogenetic core genome snp analysis, but excluded this strain from visualization of the phylogenetic core genome snp tree to get a better resolution of the cff and cfv branches. snp discovery was focused on a comparison between cff and cfv strains; therefore, cft strain 03- was replaced by cff strain 82- as a reference for this analysis. recombination regions in the core genome alignment were detected and visualized using gubbins. a phylogenetic tree was constructed using fasttree with a generalized time-reversible model and gamma correction on the recombination-filtered snps in the core genome of all isolates, including cft strain 03- the resulting tree was rooted on cft strain 03- using dendroscope prior to visualization using itol. beast analysis recombination-filtered non-synonymous snps of the mammal-associated c. fetus isolates were extracted from the gubbins results and used for divergence dating in beast, using the isolation dates as tip dates in the phylogenetic tree, as outlined in the beast manual with the following modifications:,000× sampling and a general time reversible model plus gamma correction as the distance model. combinations of a strict clock, log-normal clock, exponential clock and random local clock as the clock model and a constant population, exponentially-growing population, and a bayesian skyline plot with six groups as demographic models were used. maximum ess values and lowest 95 % confidence intervals in the predicted divergence dates of the clades were obtained with a random local clock and a bayesian skyline plot with four groups as the demographic model. calculation of branch-specific dn/ds ratios to calculate the dn/ds ratio per branch in the phylogenetic tree, we aligned the genomes of the mammal-associated c. fetus isolates with parsnp v and cft strain 03- as outgroup. synonymous and non-synonymous snps were determined on the basis of their location in the coding regions of the cff strain 82– reference genome. recombination regions were detected using gubbins and excluded from the alignment. ancestral state reconstruction of the node sequences was performed using fastml with a generalized time-reversible model and gamma correction. the dn/ds ratios were determined per branch between node sequences. comparison of snps and genes with phenotypic characteristics of the strains as traditional differentiation of c. fetus subspecies are based on the biochemical 1 % glycine tolerance and h2s production tests, the genomes were screened for genes and snps that were associated with these phenotypic characteristics, and tested for 1 % glycine tolerance and h2s production in cysteine-rich medium as described above. the protein-encoding gene presence and absence were determined using blast-based all vs all comparisons with prokka-annotated genomes using roary, which clustered the proteins using mcl-edge. a fisher’s exact test was used to calculate the two-tail probability value of respectively the detected snps and genes versus the outcome of the biochemical tests. calculation of clinical association a fisher’s exact test was used to calculate the two-tail probability value of snps and genes which were specifically found in strains that were isolated from bovine abortions. significantly-associated snps or genes were checked for their presence in cff strain 04/ as this is a non-cfv bovine clinical strain. snps and genes that were not present in this strain were excluded from this analysis, to separate clinical-associated from phylogenetically-associated snps and genes. availability of data genome sequences are available from the european nucleotide archive and from ncbi genbank, with the accession numbers listed in table   phylogenetic analysis of the core genome snps phylogenetic analysis, based on core genome snps using cft as the outgroup, showed that the mammal-associated cff and cfv genomes group into five distinct clades. the division of clades was consistent with the serotypes of the strains; clades and consist of serotype b strains and clades and consist of serotype a strains. the division of clades was also consistent with the classification of mlst sequence types, except for cff strain h1-uy, which has the cfv-associated mlst st- genotype. interestingly, snp phylogeny showed the divergence of clade and clade from a common cff ancestor. clade consists of phenotypically-identified cff, cfvi and cfv strains without a clear separation in the phylogeny. however, when the genotypic characterization is used, clade consists exclusively of cfv and cfvi strains and clade 1– of cff strains.fig. ml tree of c. fetus core genome and mammalian associated c. fetus recombination regions. cft 03- was excluded from the visualization. tree branches are labeled with the clade numbers, serotype and mlst sts of the c. fetus strains. subspecies identification, h2s production test results, the presence of the putative cysteine transporter and the gubbins predicted recombination regions in the mammalian associated c. fetus core genome are shown, with on the right side a zoomed ml tree of clade gene flge encodes a flagellar hook protein and cjei a type ii r-m system protein. scale bar at the bottom left represents the nucleotide substitutions per site beast analysis determined that diversification of the mammal-associated c. fetus is a relatively recent event under disruptive selection pressure we used beast to determine the divergence date of the different clades. the topology of the beast tree was slightly different than the maximum likelihood snps tree, especially for the serotype b strains. this may be either the result of using only synonymous snps in the beast analysis or because of the differences between ml methods and methods based on coalescence. the 95 % confidence interval of the beast analysis was too large to estimate the divergence date of clade and clades and and the divergence date of clades and. interestingly, the beast analysis showed that all currently circulating strains in clades and diverged very recently from each other.fig. beast tree with divergence dates. beast analysis 95 % ci divergence dates are shown within brackets and with blue bars on branch nodes. dn/ds ratios are shown in italics on the branches. branches are colored according the clock rates, represented as substitutions per site per year. scale bar at the bottom represents thousand years ago cfv is restricted to the bovine genital tract, but cff can be isolated from different animal species, e.g. cattle, sheep, goats, pigs, horses, and humans. therefore, we attempted to investigate whether this change in niche preference has resulted in increased evolutionary pressure on cfv. we estimated the rate of molecular evolution by comparing the ratio of non-synonymous versus synonymous snps. the dn/ds ratio was the highest in the strains of clade compared to the ratio of the strains of clade and in all branches preceding, based on node sequences reconstructed using fastml. an increased dn/ds ratio is indicative for recently diversifying genomes, whereas lower dn/ds ratios correspond to older snps and more diverged genomes. all clades have a dn/ds ratio > , with the highest dn/ds ratio of clade of, showing that the strains of clade are evolving under relatively more diversifying selection than the strains of the other clades. the estimated molecular clock rate for the mammal-associated c. fetus strains was –  × 10− substitutions per kb per year, which corresponds with the estimated molecular clock rate of c. jejuni of –  × 10− per kb/year. clade-specific snps and orthologs the numbers of snps and orthologous genes specifically present in each clade are presented in table  as well as the specific snps and orthologs for clades representing the serotypes and subspecies of the strains. the details and annotation of the detected snps and orthologs are listed in additional file 1: table s table core genome snps specific for c. fetus clades cff, c. fetus subsp. fetus; cfv, c. fetus subsp. venerealis; cfvi, c. fetus subsp. venerealis biovar intermedius clades 1– consisted only of cff serotype b strains. for clade snps and orthologs were found, specifically present in clade strains. the two clade 1-specific orthologs were hypothetical proteins with unknown function. specific for clade snps and orthologs were found. the specific orthologs were located in a region encoding mainly hypothetical proteins, as well as a fic domain protein and a transposase. this clade 2-specific region is chromosomally located in cff strain 04/. all strains of both clades and contain specific snps and specific orthologs. the clade and 2-specific orthologs included a toxin/antitoxin system, restriction-modification system-associated genes and s-layer associated genes. clades and consist of cff serotype a strains. clade contains specific snps and six specific orthologs. the specific orthologs encoded hypothetical proteins of unknown function, as well as r-m type i and r-m type iii genes. for clade specific snps and no specific orthologs are found. all strains of clades 3– contain snps and specific orthologs that are mainly located in one region that includes a glycosyltransferase, a methyltransferase and a ketoreductase. clade consists of both genotypically-identified cfv and cfvi strains. for this clade, specific snps were found. the genotypically-identified cfv strains of clade contain five specific snps and six specific orthologs. the specific orthologs encoded hypothetical proteins, a transcriptional regulator, a resolvase and a transposase. analysis of snp regions that are subject to recombination interestingly, 42 % of the snps differing between clades 1– and clade are located in regions <1 kb apart. this is suggestive of recombination, since on average a snp would be expected to occur approx. every  kb, based on the number of snps identified here in genomes of approx. mbp.. we predicted potential recombination regions in the mammal-associated c. fetus core genome alignment with gubbins and visualized the results. using this analysis, one region showed recombinations in all strains of clades 1– and three recombination regions were found in all strains of clade the recombination region common to all strains of clades 1– encoded the flagellar hook gene flge. all strains of clade contained a recombination region with a type ii r-m gene and genes involved in sugar biosynthesis, including: a nad-dependent epimerase/dehydratase; a nucleotide sugar dehydrogenase; and a polysaccharide biosynthesis protein. another clade 5-specific recombination region is cff8240_1473-cff8240_ that encodes: leader peptidase a; ribose-phosphate pyrophosphokinase; a subunit of aspartate carbamoyltransferase, a formate dehydrogenase subunit and a glutamate synthase subunit. in all strains of clade a recombination region was found, encoding: several radical sam domain proteins; two transketolase subunits; a methyltransferase; and two glycosyltransferases. additionally, the s-layer region was also identified by gubbins as a recombination region, although the sap locus itself is not included in the core genome alignment due to assembly issues of this region. deletion of a putative cysteine transporter is associated with h2s production-negative strains currently, the oie-prescribed biochemical tests to differentiate cff, cfv and cfvi consist of determining tolerance to 1 % glycine and production of hydrogen sulfide from l-cysteine. we analyzed the genomes of the c. fetus strains for the presence of specific genes associated with these phenotypic characteristics. interestingly, h2s negative strains did not encode two subunits of an amino-acid abc transporter, that together with an atp-binding subunit putatively form an amino-acid abc transport system involved in cysteine transport. the absence of the two abc transporter genes in the c. fetus h2s negative strains possibly explains the phenotypic characteristic of these strains. cfvi strain 98/ was described before as h2s positive; however, in this study, the isolates that were sequenced were biochemically characterized and the sequenced isolate of strain 98/ was h2s production negative, in accordance with the absence of the putative cysteine transporter in this genome. conjugative transfer region associated with strains from bovine abortions the c. fetus strains were arranged according to their clinical features: nine strains were from bovine abortions; seventeen strains were from screenings; and for fifteen strains, the clinical features were unknown. no snps were found that were both present in the only known cff strain from a bovine abortion and significantly present in the cfv/cfvi strains from bovine abortions. however, one region was significantly present in all strains from bovine abortions; this region contains genes encoding a filamentation induced by cyclic amp domain protein, a dna-binding protein and conjugative transfer proteins. this region is located in cff strain 04/ cfv strain 97/ and cfv strain 84- on a megaplasmid or an extra-chromosomal element. this tra/trb region is present in c. fetus strains isolated from non-abortion cases; cff strain 98/v and cfv/cfvi strains ccug 642- adri b 84- lmg nctc zaf adri and 99-. discussion the core genome snp analysis of the mammal-associated c. fetus strains identified a large number of snps. since available pcr methods and mlst typing schemes are not able to identify c. fetus strains correctly to subspecies level, cfv/cfvi-specific snps could be used to develop new diagnostic methods, such as a snp probe-pcr, or to improve the current mlst scheme. however, two drawbacks of a pcr based on cfv/cfvi-specific snps are that the stability of such snps is unknown and that the design of a reliable pcr based on specific snps requires a high degree of optimization. in a previous study, we presented the division of the mammal-associated c. fetus strains into two different clusters, based on a core gene alignment of five cff and cfv/cfvi strains. in this study, cff strains and cfv/cfvi were included, which gave a better resolution of the branches with cff strains, when compared to a previous analysis with five cff strains. the current core genome snp analysis separated the c. fetus strains into five different clades, that contain either the genotypically-identified cff or cfv/cfvi strains. we demonstrated that the biochemical differentiation of the cff, cfv and cfvi strains is not supported by the core genome snp phylogeny: clade consists of phenotypically-identified cff, cfv and cfvi strains without a congruent separation in the phylogeny of these strains this showed that the phenotypic separation of the mammal-associated c. fetus strains is not supported by the core genome snp phylogeny, which gives rise to the consideration if the current phylogenetic subspecies differentiation is still reasonable, as previously mentioned. forces of cff and cfv diversification phylogenetic analysis indicated that cfv/cfvi clade and cff clade have a common cff ancestor. no association of the clades with geographic origin and/or host specificity was observed. interestingly, the short terminal branches and recent diversification suggest that the isolates within each clade, except for strain h1-uy, have diverged recently. this is in contrast to what is observed in cft, where terminal branches of the isolates are much longer and a much higher genome diversity is observed. the world-wide spread of the mammal-associated c. fetus strains is potentially associated with the improvements in cattle breeding by selection and cross-breeding which started around the 1700s and 1800s, followed by extensive spread of high-producing dairy cows and beef cattle, which may have carried only a very limited number of c. fetus clones. we suggest that this spread from a limited number of sources has resulted in our current observation of the highly clonal nature of mammal-associated c. fetus. alternatively, the observation of only recently diverged strains is a result of selective sweeps, when a population member with an advantageous trait will take over the population before the trait can spread to other members. it is unlikely that selective sweeps occurred in the mammal-associated c. fetus population, since cattle are intensively monitored worldwide for the presence of c. fetus subsp. venerealis and only one non-clonal c. fetus strain was found. cff strain h1-uy was isolated from the blood of a rural worker in who was diagnosed with cellulitis and was in daily contact with cattle. the rates of changes at non-synonymous and synonymous snps indicate whether a gene is under purifying or diversifying selection and is expressed in the dn/ds ratio. the dn/ds ratio of bacterial genes under stabilizing selection falls generally within the range –. in the mammal-associated c. fetus, the dn/ds ratios were all > , showing that all genomes are under diversifying selection. the highest dn/ds ratio was found for the clade strains. this high ratio may be due to the genetic features of these strains: they have no crispr-cas systems and the clade 5-strains have specific snps in the type ii r-m system gene cjei compared to the other strains. type ii r-m systems can undergo major changes in specificity by recombination events as shown for helicobacter pylori and r-m type ii systems may play a role in plasmid transformation, as described for c. jejuni where knockout mutagenesis of gene cjei resulted in a strain with a -fold-enhanced transformation efficiency. it is unknown whether these snps influence the functionality or specificity of this type ii r-m system in c. fetus strains, but the clade 5-strains are possibly more susceptible to insertion of foreign dna, since the genomes contain genomic islands encoding t4sss, phages and insertion sequences. h2s negative phenotype of cfv strains associated with loss of a putative cysteine transporter the original classification of the mammal-associated c. fetus subspecies cff and cfv are based on the 1 % glycine tolerance test. cfv strains can be discriminated from cfvi strains with the h2s production test. we show that the h2s negative c. fetus strains have lost a putative cysteine transporter. without the encoded transporter, the cells are possibly less capable of importing cysteine, which under normal situations is reduced while forming h2s. the h2s-negative cfv strains have a niche restriction to the genital tract of cattle, whereas the h2s-positive cfvi and cff strains are not restricted to the genital tract and assumed to be able to colonize the intestines as well. one may speculate that the cfv strains have a defect causing the incapability to grow outside of the genital tract, but it is unknown if this partial deletion of the putative cysteine transporter can be associated with the niche restriction of these strains. virulence-associated genes specific for clades and strains from bovine abortions since the first description of the c. fetus subspecies in it is presumed that cfv cause disease in the genital tract of cows, like enzootic venereal sterility and abortions in pregnant cows, and that cff only cause sporadic abortions. it is unknown which genomic features are responsible for the pathogenicity of the c. fetus strains, but potential candidates are the surface layer of the c. fetus cells and type iv secretion systems. we studied the genomes of nine strains from bovine abortions and observed that, in addition to the s-layer proteins and t4sss, all the strains from bovine abortions contain a region encoding conjugative transfer proteins. this tra region is located in the closed genomes of strains 04/ 97/ and 84- on a plasmid/ice, which also contain a trb-t4ss and fic genes. the fic domain proteins encoded by this region are of potential interest, because fic domain proteins have an immune-modulatory function by influencing the cytoskeleton organization of the host cells and are translocated by t4sss in legionella pneumophila, coxiella burnetii and bartonella henselae, which has also been proposed for c. fetus strains. the trb/tra genes are also present in c. fetus strains isolated from non-abortion cases. it is possible that these strains can also cause bovine abortions; however, more studies, specifically animal experiments investigating the virulence of c. fetus strains are required to assess this hypothesis.  phylogenetic core genome snp analysis divided the mammal-associated c. fetus strains into five different clades, which were consistent with the serotypes, but not with the phenotypes of the strains. beast analysis showed that the clade with genotypically-identified cfv/cfvi strains has evolved from a cff ancestor under diversifying selection. phylogenetic analysis of the core genome snps did not differentiate h2s-negative cfv from h2s positive cfvi strains. the partial deletion of a putative cysteine transporter is observed in all h2s negative cfv strains. additional files additional file 1: table s table with details of c. fetus snps and genes, specific for clades and bovine abortion strains. additional file 2: figure s schematic representation of the putative cysteine transporter genes cff8240_ - cff8240_ in strain cff with flanking genes rard and psga. the start and stop positions of respectively rard and psga are shown in italics. cfv strain 97/ lacks one complete orf and contain a remnant of one orf of the putative cysteine transporter. shown are the similar deletion sites sequences in the cfv strains containing the incomplete transporter. abbreviations bgcbovine genital campylobacteriosis c. fetuscampylobacter fetus cffc. fetus fetus cftc. fetus testudinum cfvc. fetus venerealis cfvic. fetus venerealis biovar intermedius gigenomic island kyathousand years ago snpsingle nucleotide polymorphism t4sstype iv secretion system 
the prescribed methods for subspecies identification of the cff and cfv isolates are: tolerance to 1 % glycine and h2s production. campylobacter fetusbovine genital campylobacteriosissubspecies differentiationcore genome snp analysish2s productioncysteine transporterissue-copyright-statement© the author 2016 c. fetus has been divided into three subspecies: c. fetus subsp. this suggests the need for a closer evaluation of the current c. fetus subspecies differentiation, considering that the phenotypic differentiation is still applied in bgc control programs. subspecies identification of mammal-associated c. fetus strains is crucial in the control of bovine genital campylobacteriosis, a syndrome associated with cfv. multiple c. fetus clades were observed, which were not consistent with the biochemical differentiation of the strains. in this study, we observed the deletion of a putative cysteine transporter in the cfv strains, which are not able to produce h2s from l-cysteine. phylogenetic reconstruction of the core genome single nucleotide polymorphisms within cff and cfv strains divided these strains into five different clades and showed that the cfv clade and a cff clade evolved from a single cff ancestor.
 campylobacter fetus can cause disease in both humans and animals. c. fetus has been divided into three subspecies: c. fetus subsp. fetus, c. fetus subsp. venerealis and c. fetus subsp. testudinum. subspecies identification of mammal-associated c. fetus strains is crucial in the control of bovine genital campylobacteriosis, a syndrome associated with cfv. the prescribed methods for subspecies identification of the cff and cfv isolates are: tolerance to 1 % glycine and h2s production.  in this study, we observed the deletion of a putative cysteine transporter in the cfv strains, which are not able to produce h2s from l-cysteine. phylogenetic reconstruction of the core genome single nucleotide polymorphisms within cff and cfv strains divided these strains into five different clades and showed that the cfv clade and a cff clade evolved from a single cff ancestor.  multiple c. fetus clades were observed, which were not consistent with the biochemical differentiation of the strains. this suggests the need for a closer evaluation of the current c. fetus subspecies differentiation, considering that the phenotypic differentiation is still applied in bgc control programs. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords campylobacter fetusbovine genital campylobacteriosissubspecies differentiationcore genome snp analysish2s productioncysteine transporterissue-copyright-statement© the author 2016 