ABSTRACT
Ornithobacterium rhinotracheale is a causative agent of respiratory tract infections in avian hosts worldwide but is a particular problem for commercial turkey production. Little is known about the ecologic and evolutionary dynamics of O. rhinotracheale, which makes prevention and control of this pathogen a challenge. The purpose of this study was to gain insight into the genetic relationships between O. rhinotracheale populations through comparative genomics of clinical isolates from different U.S. turkey producers. O. rhinotracheale clinical isolates were collected from four major U.S. turkey producers and several independent turkey growers from the upper Midwest and Southeast, and whole-genome sequencing was performed. Genomes were compared phylogenetically using single nucleotide polymorphism (SNP)-based analysis, and then assembly and annotations were performed to identify genes encoding putative virulence factors and antimicrobial resistance determinants. A pangenome approach was also used to establish a core set of genes consistently present in O. rhinotracheale and to highlight differences in gene content between phylogenetic clades. A total of 1,457 nonrecombinant SNPs were identified from 157 O. rhinotracheale genomes, and four distinct phylogenetic clades were identified. Isolates clustered by company on the phylogenetic tree, however, and each company had isolates in multiple clades with similar collection dates, indicating that there are multiple O. rhinotracheale strains circulating within each of the companies examined. Additionally, several antimicrobial resistance proteins, putative virulence factors, and the pOR1 plasmid were associated with particular clades and multilocus sequence types, which may explain why the same strains seem to have persisted in the same turkey operations for decades.
IMPORTANCE The whole-genome approach enhances our understanding of evolutionary relationships between clinical Ornithobacterium rhinotracheale isolates from different commercial turkey producers and allows for identification of genes associated with virulence, antimicrobial resistance, or mobile genetic elements that are often excluded using traditional typing methods. Additionally, differentiating O. rhinotracheale isolates at the whole-genome level may provide insight into selection of the most appropriate autogenous vaccine strain, or groups of strains, for a given population of clinical isolates.
INTRODUCTION
Ornithobacterium rhinotracheale is a Gram-negative bacterium that causes respiratory disease in avian populations worldwide. Infections with this emerging pathogen have been reported in multiple avian host types, including both domestic poultry and wild birds (1). However, the clinical course of O. rhinotracheale is particularly severe in commercial turkeys, as it often manifests as respiratory disease complex and can lead to reductions in feed and water intake and increased flock mortality (2–4). Coinfections with O. rhinotracheale and other respiratory pathogens, such as avian influenza virus, Newcastle disease virus, avian metapneumovirus, and infectious bronchitis virus, are common and can result in higher condemnation rates than either pathogen alone (5–8).
Antibiotics are used for treatment of O. rhinotracheale in commercial turkeys, with tetracycline, erythromycin, and penicillin use common in the United States (9). However, the success of antibiotic therapy is often limited due to antibiotic resistance. A study conducted on turkeys in Minnesota demonstrated an increase in tetracycline resistance from 25% to 43% over a 6-year period (10). A more recent study in Japan revealed 100% resistance of 21 clinical isolates to amikacin, colistin, gentamicin, kanamycin, neomycin, polymyxin B, streptomycin, and sulfamethoxazole (11). However, it should be acknowledged that there is no standard protocol for antimicrobial susceptibility testing in O. rhinotracheale, making classification of isolates as susceptible or resistant highly subjective and limiting comparisons between studies. These inconsistencies emphasize the need for additional preventative measures for O. rhinotracheale that would reduce the need for antibiotic therapy.
The United States Animal Health Association recently ranked O. rhinotracheale infection as the third most economically important issue facing the turkey industry and the second most important respiratory condition after colibacillosis (12). Similarly, a survey of poultry veterinarians in Canada revealed that O. rhinotracheale is the third most frequently diagnosed bacterial infection in turkeys, behind avian pathogenic Escherichia coli and staphylococcal infections (13). Many factors have been proposed to be responsible for the increasing importance of O. rhinotracheale in turkey production, including a lack of commercially available live vaccines (14). Without commercial vaccines available, producers must rely on autogenous vaccines produced from flock-specific bacterial cultures for O. rhinotracheale prevention (15). However, autogenous vaccines in general have not demonstrated long-term effectiveness in cross-protection against different strains of the same pathogen in poultry (16, 17). Conversations with poultry veterinarians have indicated similar experiences with O. rhinotracheale (18). This may potentially be due to the limited capability of a single autogenous vaccine strain to protect the host against diverse strains of the same pathogen. Understanding the genomic differences between O. rhinotracheale strains circulating in different commercial turkey operations may provide insight into the limited efficacy of autogenous vaccines.
There are large knowledge gaps regarding the genomics of O. rhinotracheale, as few studies have examined O. rhinotracheale at the whole-genome level, and only three O. rhinotracheale genomes were publicly available at the time of this study. Most studies have used multilocus sequence typing (MLST) to classify isolates, which has comparatively low resolution and excludes the identification of genes that could be responsible for functional differences between strains (19, 20). The objective of this study was to better understand the genomic landscape of O. rhinotracheale in commercial turkey production through comparative genome sequencing of clinical isolates from different U.S. turkey producers and to use these isolates to establish an O. rhinotracheale pangenome that will distinguish acquired genetic components from those consistently present in O. rhinotracheale. The strain-to-strain differences among O. rhinotracheale isolates in this study may be used to optimize selection of the best autogenous vaccine candidate, or candidates, in order to reduce the incidence of clinical disease from currently circulating O. rhinotracheale strains.
RESULTS
This study included 157 clinical O. rhinotracheale isolates from four commercial turkey producers (designated A, B, C, and D) and two independent turkey growers. There were 27 O. rhinotracheale isolates from company A, 53 from company B, 34 from company C, 41 from company D, and two from the independent turkey growers. Collection dates for the clinical isolates ranged from 1999 to 2019, although exact collection dates were not available for all isolates. The 157 whole-genome sequences in this study had a mean sequencing depth of 55.2, a genome fraction of 88.7%, and a length of 2.3 million base pairs.
SNP-based phylogeny reveals four distinct phylogenetic clades.The single nucleotide polymorphism (SNP) analysis identified 1,457 core SNPs from 157 whole-genome sequences using the transversion model of DNA substitution with equal base frequencies. All isolates fell within four distinct phylogenetic clades (Fig. 1). These four clades were separated by a minimum of 150 SNPs and a maximum of 471 SNPs, while within these clades, individual isolates were separated by a minimum of 0 SNPs and a maximum of 72 SNPs. Clades 1 (n = 47) and 2 (n = 6) were composed exclusively of isolates of sequence type 1 (ST1), and clades 3 (n = 49) and 4 (n = 57) were composed exclusively of ST34 and ST35 isolates, respectively.
SNP-based phylogenetic tree of clinical O. rhinotracheale isolates from U.S. turkey producers. The tree represents 1,457 core SNPs and is labeled by collection date and multilocus sequence type. The different colored circles represent the different turkey producers that participated in the study. The scale bar represents 200 core SNP differences.
Clustering by production company was visually evident, and each company had isolates fall into at least two phylogenetic clades. Within each phylogenetic clade, collection dates varied widely. Clade 1 was composed of isolates ranging in collection dates from 1995 to 2018, clade 2 had isolates with collection dates of 2017 and 2018, clade 3 had collection dates ranging from 2009 to 2019, and clade 4 had isolates with collection dates from 1999 to 2019. Interestingly, within some of these clades, isolates from the same company had collection dates 20 years apart, indicating that the same O. rhinotracheale strains have persisted in these commercial turkey operations over time and have continued to be a causative agent of respiratory disease, despite various preventative efforts.
There also did not appear to be major differences by geographic region. Farms from companies A, B, and D were located primarily in the Midwest, while company C was located in the Southeast, demonstrating that geography has not been a barrier to the circulation and dissemination of these different O. rhinotracheale strains.
The core O. rhinotracheale genome contains several virulence factors and antimicrobial resistance determinants.The O. rhinotracheale pangenome consisted of 4,927 genes, of which 1,399 (28.4%) were core genes, meaning they were present in 99% of the O. rhinotracheale genomes. There were 319 (6.5%) genes present in 95 to 99% of O. rhinotracheale genomes, 943 (19.1%) present in 15 to 95% of genomes, and 2,266 (46.0%) present in 0 to 15% of genomes, demonstrating a large number of accessory genes present only in some O. rhinotracheale genomes but not others. All genes in the O. rhinotracheale pangenome are listed in Data Set S1 in the supplemental material.
Several virulence factors were found to be encoded in the core genome. This included a predicted IgM protease, which has recently been associated with IgM degradation in the early host immune response, and toxin-antitoxin biofilm protein (TabA) (21). Proteins associated with gliding motility were also identified, including gliding motility lipoprotein (GldH) and motility protein B (MotB), which serve as an important virulence mechanism in related Flavobacterium species (22, 23).
Several predicted chromosomal antimicrobial resistance proteins were also found to be encoded in the core O. rhinotracheale genome, including a macrolide export protein (MacA) and macrolide export ATP-biding/permease protein (MacB). Penicillin-binding proteins were also identified (PBP1a and PBP4), as well as multidrug resistance proteins (MdtA, MdtN, NorM, and YheL).
Gene content analysis demonstrates differences between phylogenetic clades.Several proteins were found to be unique to specific phylogenetic clades. The genomes in clade 1 contained 21 unique proteins that were absent in genomes from all other clades. Notably, these proteins included those predicted as conjugative transposons TraK, TraM, TraN, and TraO, an IS110 family transposase, and a ParB-related ThiF family protein, suggesting a potential role for mobile genetic elements in genomes from this clade, although no plasmid replicons were identified. Genomes in clade 2 had 135 unique proteins that were absent from all other clades. While the majority of these were classified hypothetical proteins, the putative virulence factor OatA was identified, which involves O-acytelation of the peptidoglycan layer, conferring resistance to antimicrobial compounds that target the cell wall, such as beta-lactam drugs, lysozymes, and bacteriocins (24). The TetR/AcrR family transcriptional regulator was also present in genomes from this clade. Clade 3 had 33 unique proteins, including a metallo-beta-lactamase, potentially conferring the ability to hydrolyze beta-lactam antibiotics (25). There were 47 unique proteins identified in clade 4, including an endoglycosidase, which allows bacteria to interfere with host glycosylation and circumvent the immune system (26). An additional protein associated with gliding motility, GlbD, was also found to be encoded in genomes of this clade. It is possible that the presence of these proteins could result in functional differences in the virulence mechanisms of different O. rhinotracheale strains. A list of unique proteins for each clade can be found in Data Set S2.
There were also multiple antimicrobial resistance proteins identified that were not unique to a specific clade, or uniformly present, but appear to be associated with a particular clade. For example, a predicted erythromycin resistance protein, ErmF, was identified in 100% (6/6) of genomes from clade 2, but outside this clade, it was only found in 5.3% (3/57) of isolates from clade 4. Similarly, a predicted macrolide resistance protein, Mef(En2), was found in 53.2% (25/47) of O. rhinotracheale genomes from clade 1 but was not identified in any other O. rhinotracheale genomes. It is therefore a possibility that O. rhinotracheale strains from different clades may respond differently to antibiotic therapy based on the presence or absence of certain antimicrobial resistance proteins.
Plasmid identification reveals strain specificity.The O. rhinotracheale plasmid pOR1 was identified in 57 of 158 (36.1%) O. rhinotracheale genomes. The 57 genomes in which the plasmid was identified all belonged to ST35, or clade 4, on the phylogenetic tree. The plasmid was not identified in any other O. rhinotracheale genomes in this study, and no other plasmid sequences were identified in any of the O. rhinotracheale genomes. The pOR1 plasmid in genomes from this study included several proteins predicted to be associated with resistance to heavy metals, including CopA, a cadmium/cobalt/zinc/H+-K+ antiporter, mercury transporter, copper-exporting ATPase, and DUF305 domain-containing protein of the ferritin protein superfamily (Fig. 2). The plasmid also contained a virulence-associated protein D (VapD) that has not been extensively characterized but has been previously associated with biofilm growth in other bacterial species (27).
Plasmid pOR1. The 17,676-bp plasmid is annotated by corresponding protein sequences, and the different colored arrows represent different protein categories.
DISCUSSION
The purpose of this study was to fill current knowledge gaps surrounding the genomics of O. rhinotracheale to better understand the ecology and evolution of this important respiratory pathogen in commercial turkey production. All 157 clinical O. rhinotracheale isolates used in this study fell into four distinct phylogenetic clades and, within these clades, demonstrated high levels of genomic similarity. This observation is consistent with previous studies that have shown low genetic heterogeneity among O. rhinotracheale isolates (1, 20). In addition to the similarities and differences seen visually, genomes from each of these clades were statistically associated with specific virulence factors or antimicrobial resistance genes that may play a role in the dissemination and persistence of these O. rhinotracheale strains in commercial turkey operations over time.
While all O. rhinotracheale genomes in each clade belonged the same multilocus sequence type, two of these clades, separated by >200 SNPs, both contained genomes belonging to ST1, indicating that this 7-gene subtyping method may not always correspond to differences seen at the whole-genome level. ST1 is the most frequent sequence type in the PubMLST database, with submissions from the United States, European countries, South Africa, and China, indicating the widespread dissemination of this sequence type, and hosts included both chickens and turkeys (28). ST34 and ST35 each only have two isolates in the PubMLST database, all from turkeys in the United States. More O. rhinotracheale genomes from multiple host types and geographic regions are needed to fully explore the relationship between host, geography, sequence type, and accessory gene content.
Genomes from each phylogenetic clade were associated with specific virulence and antimicrobial resistance genes that may provide insight into different ways this fastidious bacterium can successfully colonize the respiratory tract, evade the host immune response, and persist after the application of antimicrobial drugs. Eventually, a clade-specific PCR assay may be able to provide crucial information to producers, veterinarians, and researchers on the potential success of different antimicrobial compounds when access to whole-genome sequencing (WGS) is limited. For example, the O. rhinotracheale genomes in clade 2 encoded TetR and EmrF proteins, conferring resistance to tetracyclines and erythromycin, respectively, as well as OatA, associated with beta-lactam resistance. Therefore, identification of a clinical isolate belonging to clade 2 in a commercial turkey barn may lead to the application of alternative antimicrobials that have a larger potential for success in the treatment of O. rhinotracheale.
The pOR1 plasmid, identified in all genomes belonging to clade 4, may contribute to antimicrobial resistance of these O. rhinotracheale isolates through mobile genes conferring resistance to heavy metal compounds. The potential for coselection of antimicrobial resistance along with tolerance to heavy metals has been previously described (29, 30). In particular, zinc and cadmium are the most commonly observed heavy metals associated with antimicrobial resistance, and the pOR1 plasmid contained a cadmium, cobalt, and zinc antiporter, among other heavy-metal-associated proteins (31). It is possible that this plasmid may play a role in coselection of resistance to antibiotics often used to treat O. rhinotracheale.
The main limitation to this study, however, is that phenotypic data on antimicrobial resistance for these strains were not available. Therefore, the identified antimicrobial resistance elements represent a potential for resistance, rather than an exact reflection of antimicrobial resistance seen in the field. However, phenotypic observations were beyond the scope of this genomics-based study, and the lack of specific cutoff values for susceptibility or resistance of O. rhinotracheale isolates would have prevented accurate and objective classification of isolates had this been attempted. Additional investigation is needed in order to develop MIC breakpoints for O. rhinotracheale and determine the similarities and differences between genotypes and phenotypes.
This study utilized comparative genomics of O. rhinotracheale at the whole-genome level and included O. rhinotracheale isolates from some of the largest commercial turkey producers in the United States. The identification of O. rhinotracheale genomes within the same company belonging to different phylogenetic clades indicates that each commercial turkey company has multiple circulating strains responsible for clinical disease, which has implications for selection of autogenous vaccine strains. Therefore, it may be necessary to include multiple O. rhinotracheale isolates in an autogenous vaccine, or to rotate the use of individual isolates, for complete cross-protection against different O. rhinotracheale strains. The identification of genomically similar O. rhinotracheale isolates across multiple turkey producers may be due to the vertical integration of the poultry industry. Although this study did not include O. rhinotracheale isolates from matched breeder flocks, vertical transmission of O. rhinotracheale has been previously demonstrated (32). Further studies are needed to fully assess the relationship between clinical isolates and autogenous vaccines, as well as potential transmission dynamics, but this study serves as an initial step toward a better understanding of O. rhinotracheale genomics and the distribution of different O. rhinotracheale strains in commercial turkey operations in the United States.
MATERIALS AND METHODS
This was a retrospective cross-sectional study using O. rhinotracheale isolates collected from U.S. commercial turkey operations from 1999 to 2019 that were submitted to the University of Minnesota Mid-Central Research and Outreach Center (MCROC, Willmar, MN) as a part of routine poultry pathogen surveillance.
Isolation of O. rhinotracheale.All O. rhinotracheale isolates examined in this study were recovered from clinical swabs collected from respiratory lesions of commercial turkeys displaying typical signs of O. rhinotracheale as identified by field veterinarians. Swabs were streaked onto 5% sheep blood agar (SBA) and incubated in a microaerophilic environment for 48 h at 37°C. Isolated colonies were selected and grown in 3 ml brain heart infusion (BHI) broth for 48 h at 37°C. Isolates were then stored at –80°C in BHI broth containing 20% glycerol. Prior to whole-genome sequencing, an O. rhinotracheale-specific PCR was utilized to confirm the isolation of O. rhinotracheale, as previously described (33).
DNA extraction and sequencing.DNA from O. rhinotracheale isolates was extracted using the Qiagen DNeasy blood and tissue kit (Valencia, CA) according to the manufacturer’s instructions. Genomic DNA libraries were created using a Nextera XT library prep kit and Nextera XT index kit v. 2 (Illumina, San Diego, CA) according to the manufacturer’s instructions. Sequencing was performed on an Illumina MiSeq system using a 2 × 250-bp dual-index approach at the MCROC Laboratory in Willmar, Minnesota.
Single nucleotide polymorphism-based analysis.The raw sequencing reads were quality checked using FastQC v. 0.11.17, and the sequencing depth was calculated using BEDTools v. 2.29.0 (34, 35). Trimmomatic v. 0.33 was used to cut Nextera adapter sequences, remove nucleotide sequences with an average quality per base of less than 15, and discard reads smaller than 36 bases long (36). Core single nucleotide polymorphisms (SNPs) were identified from the trimmed, paired-end reads against reference genome ORT-UMN 88 (GenBank accession no. CP006828.1) using Snippy v. 4.1 with default parameters (37, 38). The core SNP alignment was run through Gubbins v. 2.3.4 to remove SNPs found in recombinant regions from the alignment (39). A SNP-based maximum-likelihood phylogenetic tree was built using IQ-TREE with the ModelFinder tool and ultrafast bootstrap approximation with 1,000 iterations (40–42). The resulting phylogenetic tree was annotated using iTOL v. 5.2, and SNP differences within and between phylogenetic clades were identified with snp-dists v. 0.7 (43, 44).
Genome assembly.O. rhinotracheale genomes were assembled de novo with SPAdes v. 3.10 using default parameters (45). Contigs with less than 500 bp were removed. Genome assemblies were assessed for quality with QUAST v. 5.0.2 (46).
MLST.All genomes were sequence-typed in silico using the publicly available MLST database with the O. rhinotracheale scheme developed by Thieme et al. in order to visualize the relationship between sequence types and whole genomes on the phylogenetic tree (20, 47).
Pangenome analysis.Open reading frames for each genome were defined using Prokka v. 1.13, and Roary v. 3.10.2 was used to generate a core genome alignment and build the O. rhinotracheale pangenome (48, 49). O. rhinotracheale genomes were included in the pangenome analysis if they had a genome fraction of at least 80%, meaning that the genome assemblies for those isolates covered 80% or more of the ORT-UMN 88 genome. All isolates that did not meet this requirement were excluded because they would have decreased the number of “core” genes in the pangenome analysis. To be considered a core gene, that particular gene had to be present in at least 99% of isolates, so it was important that the included isolates cover the majority of the O. rhinotracheale genome in order to capture the most accurate representation of core genes. However, the ORT-UMN 88 genome was not included in the pangenome analysis because that would have excluded the identification of novel genes potentially present in the clinical isolates with more recent collection dates.
Gene content analysis.Scoary v. 1.6.16 was used to analyze differences in gene content between phylogenetic clades (50). Genes that were unique to specific phylogenetic clades were further characterized by searching the corresponding amino acid sequences against the BLAST protein database with a minimum percent identity of 90% and coverage of 50% for identification (51).
Plasmid detection.There is one known plasmid associated with O. rhinotracheale, pOR1 (NC_011414.1), first characterized in 2004 (52). BLAST was used to search each assembly for the pOR1 reference sequence (51). Contigs with greater than 90% identity to the pOR1 reference with at least 50% coverage were considered to contain the plasmid. Genomes were also searched against the PlasmidFinder database to determine if any known plasmid replicons were present (39, 53). Plasmid visualization was done with SnapGene Viewer (GSL Biotech).
Data availability.All raw data generated in this project are available in the NCBI Sequence Read Archive database under BioProject PRJNA524749.
ACKNOWLEDGMENTS
This study was funded through grants from the Minnesota Turkey Research and Promotion Council and the University of Minnesota MnDRIVE initiative and through USDA-NIFA Agriculture and Food Research Initiative competitive grant no. 2018-68003-27464.
Bioinformatics were supported using tools available from the Minnesota Supercomputing Institute.
FOOTNOTES
- Received 9 December 2019.
- Accepted 25 March 2020.
- Accepted manuscript posted online 3 April 2020.
Supplemental material is available online only.
- Copyright © 2020 American Society for Microbiology.
