Animal Sampling

We retrieved G. anatis isolates during a 2-year period (2017–2018) from 10 calves from 7 unrelated farms in Belgium; all 10 calves had a history of respiratory problems (≈5% of the total amount of samples). No poultry was present at these farms; however, at farm 2 (Table 1), raw eggs were occasionally fed to the calves. We obtained all isolates from animals 4–60 days old (Table 1) exhibiting signs of infectious bronchopneumonia, such as fever (>39.3°C), cough, nasal discharge, depression, and adventitious lung sounds. Before the sampling, each calf had already been treated unsuccessfully with first- or second-line antimicrobial drugs. Thoracic ultrasound examination, performed with a 7.5-MHz linear probe as described previously (16), showed a consolidated zone in the lung of >1 cm³ in all animals. A nonendoscopic bronchoalveolar lavage (nBAL) was conducted in all cases, as described previously (17). The sampling method was approved by the ethics committee of the Faculty of Veterinary Medicine, Ghent University (approval no. EC 2016/20).

Identification

We inoculated all nBAL samples on an Oxoid Columbia blood agar enriched with 5% sheep blood (http://www.oxoid.com) and on a BD Difco modified pleuropneumonia-like organism agar plate (https://www.bd.com) containing 832,000 IU/L polymyxin, 0.36 g/L ampicillin, 23.1% deactivated horse serum, and 6.5% yeast extract for the isolation of Mycoplasma spp. We incubated blood agar plates overnight and pleuropneumonia-like organism agars for 5 days, both at 35°C and in a 5% CO₂ enriched atmosphere. We identified bacterial colonies, grown on both agars, with matrix-assisted laser desorption/ionization time-of-flight mass spectrometry by using the direct transfer method and α‐cyano‐4‐hydroxycinnamic acid as matrix, according to the manufacturer’s guidelines. We considered identifications with a log score value >2.0 to be reliable at the species level. We subcultured G. anatis isolates on Columbia blood agar enriched with 5% sheep blood (Oxoid) to obtain a pure culture, which we stored at −80°C for further analysis.

Antimicrobial-Susceptibility Testing

For susceptibility testing, we performed the broth microdilution technique for ampicillin, ceftiofur, doxycycline, enrofloxacin, florfenicol, gentamicin, kanamycin, penicillin, spectinomycin, tetracycline, tilmicosin, trimethoprim/sulfamethoxazole, tulathromycin, and tylosin, according to Clinical and Laboratory Standards Institute standards (18,19). Concentrations of all antimicrobial drugs ranged from <0.03 to >128 μg/mL. We performed susceptibility testing of amoxicillin/clavulanic acid by using the gradient strip test. We used Escherichia coli ATCC 25922 and Staphylococcus aureus ATCC 29213 as quality-control strains. In addition, we included E. coli ATCC 35218 as the quality-control strain for amoxicillin/clavulanic acid testing. We used ampicillin, tetracycline, enrofloxacin, tylosin, florfenicol, spectinomycin, and trimethoprim/sulfamethoxazole as class representatives of the penicillins, tetracyclines, fluoroquinolones, macrolides, phenicols, aminocyclitol/aminoglycosides, and potentiated sulphonamides, respectively, to determine phenotypic resistance for these classes, using Clinical and Laboratory Standards Institute breakpoints for G. anatis (Appendix Table 1) (18).

Whole-Genome Sequencing

We prepared genomic DNA by using the Bioline Isolate II Genomic DNA kit (Meridian Bioscience, https://www.meridianbioscience.com), following the manufacturer’s instructions. We constructed sequencing libraries by using the Illumina Nextera XT DNA sample preparation kit and then sequenced isolates using the MiSeq Reagent v3 kit with a 250-bp paired-end protocol (Illumina, https://www.illumina.com) according to the manufacturer’s instructions. We have deposited all generated WGS data in the National Center for Biotechnology Information Sequence Read Archive (20) under accession number PRJNA541488. We cleaned and assembled raw reads (Appendix Table 2) and used Kraken 0.10.5 (21) to perform k-mer–based classification of cleaned reads against an in-house dump of the complete genomes from the National Center for Biotechnology Information RefSeq Microbial Genomes Database (22). We analyzed paired-end reads and orphaned reads (i.e., reads where only 1 read of the pair survived cleaning) separately by using default settings and then combining the results by concatenating the output files.

Antimicrobial-Resistance Genotyping

We performed genotypic resistance gene detection, as described by Bogaerts et al. (23), against the ResFinder database (24). We defined AMR gene clusters as resistance genes on the same contig within a sample. We performed detection of mutations linked with increased fluoroquinolone MICs in the quinolone-resistance determining regions of gyrA and parC by aligning these regions in the E. coli K12 reference genome in NCBI (accession no. NC_000913.3) for gyrA (accession no. NP_416734) and parC (accession no. NP_417491.1) by using the Needle tool for pairwise sequence alignment of the EMBOSS suite (https://www.ebi.ac.uk/tools/psa) (25). We used mlplasmids 1.0.0 (https://sarredondo.shinyapps.io/mlplasmids) to predict whether assembled contigs were either plasmid- or chromosome-derived, by using E. coli as species model and 1,000 bp as the minimum sequence length (26). We then compared contigs predicted to be plasmid-encoded by using blastn (https://blast.ncbi.nlm.nih.gov/Blast.cgi), with default settings, against the nucleotide database. We performed transposase detection by using ISFinder (https://www-is.biotoul.fr/index.php) with the blastn tool, using default settings (27), to substantiate the presence of transposable elements in close proximity to the AMR gene clusters in the specific contigs of the whole assembly. Last, we used ICEberg 2.0 (http://db-mml.sjtu.edu.cn/ICEberg), with default settings, to detect integrative and conjugative elements (ICEs) or integrative and mobilizable elements (IMEs) in the G. anatis assemblies (28).

Sample Relatedness

For multilocus sequence typing (MLST), we used an in-house copy of the MLST database for G. anatis hosted by the PubMLST platform (http://pubMLST.org/anatis) (29), which we pulled in-house using the REST API (30), for MLST genotyping. We typed individual loci separately by aligning the assembly for each sample against all allele sequences of that locus by using nucleotide BLAST+ 2.6.0, with default values (31). We then performed filtering and best hit identification, as described previously, for AMR gene characterization. Because MLST offered limited resolution in the relationship between samples, we used a single-nucleotide polymorphism (SNP) genotyping approach based on an in-house implementation of the CSI Phylogeny workflow (https://omictools.com/csi-phylogeny-tool) (Appendix Table 3) (32), using the NCBI RefSeq entry for G. anatis (accession no. NC_015460) as reference to compare diversity among samples. We used MEGA-Computing Core 10.0.4 (https://www.megasoftware.net) to detect the best evolutionary model and construct a maximum-likelihood phylogenetic tree on the basis of the SNP matrix, setting the following options: “missing-data” set to “partial_deletion,” “site-cov-cutoff” set to 50, “branch-swap” set to “very_weak,” “ml-method” set to “spr3,” “action” set to “model,” and “bootstraps” set to 100. We then repeated the same workflow by using the genome assembly of isolate GB8 (Appendix Table 3), filtered on contigs >1,000 bases with a k-mer coverage of 10–50× as reference. We visualized the resulting phylogenetic trees by using iTOL (33) and, afterward, a midpoint rooting. In addition, we constructed a core genome MLST (cgMLST) scheme to investigate the relationship of the isolates in Belgium compared with all genomes for this species publicly available in the NCBI database (Appendix Table 4).

Identification

We compiled all strain origin information and co-infection data (Table 1). The G. anatis isolates were all nonhemolytic and were recovered as a pure culture (50% of cases) or the predominant isolate in large numbers (50% of cases). When a dominant culture was obtained, other pathogens were detected to a lesser extent. All calves recovered from the pneumonia because of appropriate antimicrobial therapy, except 1 who was euthanized because of cardiac failure.

Antimicrobial Susceptibility Testing

We observed high MIC values for tylosin, tetracycline, spectinomycin, kanamycin, and enrofloxacin for all isolates, which most likely explains therapeutic failure (Table 2; Appendix Table 1). All isolates exhibited very low MIC values for ceftiofur and amoxicillin/clavulanic acid.

Whole-Genome Sequencing

The number of raw paired-end reads, genome assembly length, N50 (a metric used as a proxy for assembly quality that was defined as the length at which contigs of equal or longer length contained >50% of the assembled sequence), and number of contigs >1,000 bases was in the same range for all samples, with a median of 372,623 raw paired-end reads, median assembly length of 2,483,037 bases, median N50 value of 105,124 bases, and median of 58 contigs >1,000 bases across all samples (Appendix Table 2). Genome assembly sizes were close to the expected size of ≈2.69 Mb (34), indicating high quality of the WGS run. K-mer–based classification of read content for all isolates confirmed the samples to be G. anatis, given that this was the only species identified in the sample having a 5% read cutoff.

AMR Genotyping

By using the ResFinder database, we detected various AMR determinants in the WGS data for all isolates (Table 2). In total, we detected 24 different resistance genes across all 10 isolates, and several genes were present in multiple isolates. We found all isolates harbored resistance genes targeting aminoglycosides, phenicols, macrolides, sulphonamides, and tetracyclines. Seven isolates also harbored resistance genes such as bla_CARB-8 or bla_TEM-2 targeting β-lactamase–susceptible penicillins. Six isolates contained dfrA1, conferring resistance against trimethoprim. Isolate GB10 carried qnrD1, a plasmid-mediated quinolone resistance determinant. We found mutations linked with increased fluoroquinolone MICs in the quinolone resistance determining region of gyrA and parC (35) in all isolates, including a single-point mutation in parC (Ser-80 to Ile) and 2 mutations in gyrA resulting in S-83 to Y or F, and D-87 to A or G, changes. We determined the genotype to phenotype correspondence to be 90% (phenotypic observations might be explained by genotypic detection of corresponding resistance genes). In GB10, we found very high MIC values for penicillin/ampicillin and no corresponding resistance gene. We did find resistance genes without corresponding high MIC values for potentiated sulphonamides in isolate GB3 and for phenicols in isolates GB5, GB7, GB8, GB9, and GB11.

Some resistance genes were closely linked into resistance gene cassettes (Table 3). Overall, we observed a high diversity of resistance genes, both in determinants present in resistance gene clusters and in separate contigs. We detected gene clusters with 3–4 of the same resistance genes found in GB4, GB9, and GB11, and 2 identical resistance genes in GB3 and GB6 (Table 3). In 19 of 20 clusters, we observed a link with transposases in close proximity or localization on putative predicted IMEs, plasmids, or both (Table 3). In addition, we detected a type 4 secretion system not associated with a resistance gene cluster in GB2, GB5, GB7, and GB10 (data not shown).

Sample Relatedness

To evaluate the relationship between isolates, we performed MLST by using the public G. anatis database hosted by the PubMLST platform. However, an exact allelic match could only be identified for 1 locus in GB2, 2 loci in GB3, 2 loci in GB4, 3 loci in GB 5, 2 loci in GB6, 1 locus in GB7, 1 locus in GB8, 2 loci in GB9, 1 locus in GB10, and 2 loci in GB11 (in a total of 8 loci in the scheme). Reliable allele calling for the remaining loci was not possible because of mismatches and different lengths for all samples. Closer inspection revealed that the MLST database only contained 89 isolates corresponding with 81 profiles, suggesting that MLST failed because of the lack of an available background to compare against.

Because MLST was not appropriate for delineating relationships, we performed SNP genotyping by using the NCBI RefSeq reference for G. anatis (UMN179). We found 14,583–15,234 SNPs for all samples (Appendix Table 3), resulting in a total SNP matrix of 32,104 positions, indicating large diversity between samples. We repeated the workflow by using the assembly of GB8 (which had the highest original read mapping rate) as a reference; this step ensured that the number of SNPs was not erroneously inflated by taking a reference not suited for SNP genotyping (i.e., a reference too divergent from the actual samples). We found 8,978–11,137 SNPs for all samples (Appendix Table 3), resulting in a total SNP matrix of 25,166 positions, confirming the large genetic diversity among samples. Afterward, we performed model selection and phylogenetic tree reconstruction with MEGA, identifying the general time reversible model as the best fit for both references.

Figure 1

Figure 1. Phylogeny of Gallibacterium anatis isolates from cattle in Belgium, 2017–2018, based on single-nucleotide polymorphism genotyping when using GB8 as a reference. Node labels indicate bootstrap support values (expressed as decimals). Branch...

We used GB8 as reference for 1 phylogentic tree (Figure 1) and G. anatis UMN179 as reference for another (Appendix Figure). Although branch lengths differed, their underlying topology was identical and well supported by high bootstrap values, indicating that, although some isolates clustered together with fewer differences (GB10 with GB2, GB4 with GB9 and GB11, GB3 with GB6), overall we observed large variation between the different isolates. Notably, for the 4 isolates GB3, GB4, GB5, and GB6 obtained from the same farm (Table 2), only GB3 and GB6 clustered together, whereas GB4 and GB5 were located elsewhere in the phylogeny.

Figure 2

Figure 2. Phylogeny of Gallibacterium anatis isolates from cattle in Belgium, 2017–2018, based on a core genome multilocus sequence typing scheme constructed by using the 10 cattle isolates and 27 poultry isolates from...

We also constructed a cgMLST scheme on the basis of our mining all publicly available G. anatis genomes from NCBI, including in total 27 isolates from poultry, complemented with the strains from Belgium (Figure 2). Despite the existence of generally very large distances between all samples, the resulting topology indicated that the strains isolated from cattle in Belgium clustered together and were distinctly separated from all other strains isolated from poultry. Moreover, the subtopology of the isolates from Belgium was concordant with results from the SNP analysis.