Figures
Abstract
Pasteurella multocida is the primary causative agent of a range of economically important diseases in animals, including haemorrhagic septicaemia (HS), a rapidly fatal disease of ungulates. There is limited information available on the diversity of P. multocida strains that cause HS. Therefore, we determined draft genome sequences of ten disease-causing isolates and two vaccine strains and compared these genomes using a range of bioinformatic analyses. The draft genomes of the 12 HS strains were between 2,298,035 and 2,410,300 bp in length. Comparison of these genomes with the North American HS strain, M1404, and other available P. multocida genomes (Pm70, 3480, 36950 and HN06) identified a core set of 1,824 genes. A set of 96 genes was present in all HS isolates and vaccine strains examined in this study, but absent from Pm70, 3480, 36950 and HN06. Moreover, 59 genes were shared only by the Asian B:2 strains. In two Pakistani isolates, genes with high similarity to genes in the integrative and conjugative element, ICEPmu1 from strain 36950 were identified along with a range of other antimicrobial resistance genes. Phylogenetic analysis indicated that the HS strains formed clades based on their country of isolation. Future analysis of the 96 genes unique to the HS isolates will aid the identification of HS-specific virulence attributes and facilitate the development of disease-specific diagnostic tests.
Citation: Moustafa AM, Seemann T, Gladman S, Adler B, Harper M, Boyce JD, et al. (2015) Comparative Genomic Analysis of Asian Haemorrhagic Septicaemia-Associated Strains of Pasteurella multocida Identifies More than 90 Haemorrhagic Septicaemia-Specific Genes. PLoS ONE 10(7): e0130296. https://doi.org/10.1371/journal.pone.0130296
Editor: Jingfa Xiao, Beijing Institute of Genomics, CHINA
Received: February 1, 2015; Accepted: May 19, 2015; Published: July 7, 2015
Copyright: © 2015 Moustafa et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited
Data Availability: All genomic data files are available from the NCBI database; Accession numbers: JQEA00000000, JQAO00000000, JQEB00000000, JPHI00000000, JQAB00000000, JQAC00000000, JQAD00000000, JQAE00000000, JQAF00000000, JQAG00000000, JQAH00000000, JQAI00000000.
Funding: Funding was provided by the Australian Research Council Centre of Excellence in Structural and Functional Microbial Genomics, Grant number CE0562063; Australian Research Council, http://www.arc.gov.au/. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Pasteurella multocida is a Gram-negative, nonmotile, nonspore-forming coccobacillus. It is the causative agent of a spectrum of economically important diseases worldwide, including atrophic rhinitis in pigs, haemorrhagic septicaemia (HS) in cattle and buffalo, fowl cholera in poultry, snuffles in rabbits and sporadic human infections that often follow dog or cat bites [1,2]. P. multocida is a heterogeneous species with strains being commonly differentiated by serology [3], or more recently capsular locus-specific multiplex PCR [4], into five capsular serogroups designated A, B, D, E and F. Strains belonging to capsular serogroups A, D and F produce capsules composed of hyaluronic acid, heparin and chondroitin respectively [4]. The composition of the B and E capsules is unknown but the genes required for their biosynthesis have been defined [4]. A second serological typing system is also often used to differentiate strains into 16 Heddleston serotypes or serovars based on lipopolysaccharide (LPS) antigens [5]. Full strain designations usually combine both systems, such that a designation of B:2 indicates capsule serogroup B and LPS serovar 2.
Haemorrhagic septicaemia (HS) is an acute and generally fatal disease which occurs mainly in cattle and buffalo [2]. Haemorrhagic septicaemia is prevalent in Asia and Africa, where its presence is of great economic importance. In Pakistan, it has been reported as the most economically important bacterial disease of cattle and buffalo [6]. Similarly, in Thailand, HS ranks high on the list of economically important diseases of livestock [7]. Haemorrhagic septicaemia is caused by infection with P. multocida strains belonging to capsular serogroups B and E [2]. Haemorrhagic septicaemia strains producing a serogroup B capsule are predominant in Asia while strains producing a serogroup E capsule are predominant in Africa, although African serogroup B isolates and Asian serogroup E isolates have occasionally been reported [8]. P. multocida strains that cause HS belong to the LPS serovars 2 or 5 which share the same LPS outer core biosynthesis locus and produce structurally highly related, but antigenically distinct, LPS molecules [9].
The first complete genome sequence of a P. multocida strain (Pm70; GenBank accession AE004439) was determined in 2001 [10]. Analysis of the Pm70 genome identified more than 100 genes predicted to be involved in virulence and identified complete gene sets for the following pathways: TCA cycle, glycolysis, glyconeogenesis, oxidative pentose phosphate and Entner—Doudoroff [10].
There are currently 25 publicly available complete or draft P. multocida genomes. These genomes are from strains isolated from different hosts and which cause a range of different diseases [11]. There have been only limited comparative analyses of these different genomes. However, analysis of the Pm70, 36950, 3480, HN06, X73, and P1059 genomes identified a unique 18 kbp region in the porcine atrophic rhinitis isolate, HN06, that contained 14 genes, including the toxA gene encoding the P. multocida toxin (PMT) (which causes the signs of atrophic rhinitis) as well as several phage-related genes [12]. Further analyses also showed that an integrative conjugative element (ICE), ICEPmu1, was found in bovine respiratory disease isolate 36950, but not in any of the other strains. This element carried 11 different antibiotic resistance genes. A similar ICE has also been found in Histophilus somni and Mannheimia haemolytica which are both bovine respiratory pathogens [12].
Previous analysis of five P. multocida genomes (M1404, Pm70, 36950, X73, and P903) identified a core set of 1786 genes (88% of Pm70 gene content) common to all strains and a pan genome of more than 2,800 genes. Furthermore, each of these strains contained between 90 and 261 unique genes not found in any of the other strains examined [13]. For strain 36950, more than 47% of the unique genes identified were within the ICEPmu1 element, whereas for strain M1404 28% of unique genes were phage-derived elements. Importantly, a previous phylogenetic comparison of nine P. multocida strains indicated little correlation between phylogenetic relationship and disease type, capsular/LPS type, host predilection or place of isolation [13].
We recently analyzed the genotypes of 23 P. multocida isolates, 14 recovered from HS-diseased cattle and buffalo located in different geographical areas and climate zones of Pakistan and nine from different regions of Thailand. All isolates were serovar B:2 and indistinguishable by multi locus sequence typing (MLST) with all strains sequence type 122. Furthermore, all isolates from within each country were indistinguishable by pulsed-field gel electrophoresis (PFGE) [14]. Therefore, to determine whether there was any diversity across these isolates we determined whole draft genome sequences of a selection of 12 of the strains using next generation sequencing (NGS). The draft genomes were then compared with the M1404 genome (a North American serovar B:2 HS-associated strain) and the genomes from four strains not associated with HS (Pm70, 36950, 3480 and HN06). To our knowledge, this is the first detailed genomic analysis of multiple HS-associated isolates of P. multocida.
Materials and Methods
Ten P. multocida strains that had been collected from buffalo or cattle with HS, each from different regions of Pakistan or Thailand, were sourced from the National Veterinary Laboratory in Islamabad or the Department of Livestock Development in Thailand, respectively. In addition, two vaccine strains from Pakistan were also included in this study (Table 1). Except for the Faisalabad isolate (Table 1), all of the isolates had previously been identified as P. multocida and typed using MLST (all sequence type 122) [14]. The Faisalabad isolate was received from the National Veterinary Laboratory, Islamabad, Pakistan, and identified as a HS-associated P. multocida strain using both a P. multocida-specific and a HS-specific PCR (data not shown) [15].
Ethics statement
All strains were sourced from the National Veterinary Laboratory in Islamabad or the Department of Livestock Development in Thailand so no ethics approvals were required.
Whole genome sequencing
Genomic DNA was purified from each of the twelve P. multocida strains using the Qiagen DNeasy blood and tissue kit (Qiagen Cat# 69504) using 5 mL of overnight cultures grown at 37°C in brain heart infusion (BHI) broth (Oxoid, UK) and following the Gram-negative bacterial protocol outlined in the manufacturer’s instructions. DNA quantification and purity analysis was performed by agarose gel electrophoresis and Qubit Fluorometry (Life Technologies, USA). The purified genomic DNA was sequenced using the paired-end 90-bp sequencing protocol on an Illumina HiSeq 2000 (Illumina Inc., San Diego, USA) at the Beijing Genomics Institute (BGI), China. The raw read sequences were filtered to eliminate low quality reads using the following criteria; all reads with > 40% low quality (Q20) bases (parameter setting at 36 bp), > 10% Ns (parameter setting at 9 bp) or > 15 bp overlap with Illumina TruSeq adapter sequences (parameter setting at 15 bp) were removed.
Sequence assembly
The genomes of ten of the strains were de novo assembled using SPAdes v2.5.0 [16] and the remaining two strains were assembled using Velvet v1.2.07 [17]. For ten of the twelve strains, the SPAdes procedure generated acceptable assemblies. However, for the Faisal and ATTK strains, SPAdes produced unexpectedly large genome assemblies (8,278,703 and 6,134,337 bp for Faisal and ATTK, respectively) due to low level contamination with genomic DNA from another bacterial species (Bacillus cereus and Bacillus subtilis, respectively). Examination of the Velvet statistics showed that the contaminating sequences were represented by short contigs having sequencing depth below 10x. To filter this contamination, we assembled the Faisal and ATTK sequences using Velvet with a manual setting of 10 for the minimum required coverage (“velvetg-cov_cutoff 10”) to remove these undesirable components of the assembly graph prior to repeat resolution and contig extraction. Final genome assemblies of these strains were also checked for contaminating sequences using BLASTN v2.2.26 [18,19]. To evaluate the accuracy of the generated contigs of the 12 assembled genomes, they were compared with the Pm70 reference genome [10] using QUAST v2.3 [20]; sequence and assembly statistics of the 12 genomes are given in Table 2. For all strains, scaffolds (or contigs for the Velvet assembled genomes) of less than 200 bp in length were removed before the final reordering using Mauve [21] with Pm70 used as the reference sequence. The final ordered and oriented scaffold sequences were then annotated using the NCBI Prokaryotic Genome Annotation Pipeline [22]. The 12 annotated genomes were submitted to GenBank [23,24]; accession numbers are given in Table 3.
Variant calls: single-nucleotide polymorphisms (SNPs) and insertion/deletion polymorphisms (indels)
Snippy (available at https://github.com/tseemann/snippy) was used to identify SNPs and Indels in the NGS sequence reads (FASTQ format) from each of the 12 genomes compared with the reference genome, Pm70.
Genome alignments and feature analysis
Mauve v2.3.1 (default settings) [25] was used to align the genomes of the twelve Pakistani and Thai P. multocida strains with the genomes of the following strains; M1404 (bovine HS isolate, type B:2) [13], Pm70 (avian fowl cholera isolate) [10], HN06 (porcine atrophic rhinitis isolate) [26] and 3480 (GenBank: NC_017764.1) and 36950 [27] (porcine and bovine respiratory disease isolate respectively). The “homologs table” in Mauve was used to identify colinear orthologs across each of the genomes using 50% DNA identity and 50% gene coverage as the minimum criteria for a match. All genes present only in a single strain were then manually checked by BLAST [18,19] to confirm that they were unique. Genes unique to the Pakistani or Thai strains, or unique to M1404, were also identified by this method.
The PHAge Search Tool (PHAST) [28] was used to identify the positions of putative phage elements in all genomes. For PHAST analysis the FASTA files containing all concatenated contigs for each genome were uploaded to the public PHAST web server (http://phast.wishartlab.com/). To avoid false positives caused by non phage-related mobile genetic elements, PHAST filters out these mobile genetic elements using a two-step process. In the first step, PHAST uses a customized mobile genetic element database to directly filter out some of the most typical mobile genetic elements (Y. Zhou, personal communication). In the second step, PHAST filters out the rest of the mobile genetic elements when it identifies potential prophages by the density-based spatial clustering of applications with noise (DBSCAN) algorithm [28]. Specifically, the DBSCAN algorithm marks out the random mobile genetic elements as noise and clusters other gene elements into potential prophages (Y. Zhou, personal communication). In addition, PHAST predicts potential prophages based on a number of factors including the relative density of identified prophage-like genes, GC ratio, functional completeness and gene similarity to already known phages. The genomes were also checked for the presence of antimicrobial resistance genes using the ResFinder tool [29].
The capsule biosynthesis and LPS outer core biosynthetic loci were identified by BLAST comparison [18,19] against the previously reported M1404 gene clusters for capsule (20,418 bp in size) and LPS outer core (4,887 bp in size) biosynthesis. The presence and integrity of the two LPS heptosyltransferase genes, hptA and hptB were checked in order to predict whether the two different, simultaneously expressed, LPS inner core structures (glycoform A and glycoform B) were present [30]. In addition, the twelve strains were classified into either Heddleston serovar 2 or 5 by analysis of the lpt-3 gene, required for the addition of phosphoethanolamine (PEtn) to the 3 position of the second heptose (Hep II) [9].
Phylogenetic trees
The phylogenetic relationship between the strains was predicted by analysis of core genome single nucleotide polymorphisms (CG-SNPs). Identification of CG-SNPs and phylogenetic analysis was assessed using Wombac v1.2 (https://github.com/tseemann/wombac). The four closely related species, Gallibacterium anatis, Mannheimia haemolytica, P. bettyae and P. dagmatis were used as outgroups. The P. multocida strains used in these analyses were the known HS-associated strains VTCCBAA264 (GenBank: ALYC00000000) [31], P52VAC (GenBank: ALBZ00000000), Anand1C (GenBank: ALBY00000000), Anand1B (GenBank: ALBX00000000), M1404 [13] and the 12 HS isolates sequenced in this study (Table 1); the fowl cholera isolates Pm70 (GenBank: NC_002663.1), X73, (GenBank: AMBP00000000) [32], VP161 [13], P1059 (GenBank: AMBQ00000000) [32] and Anand1P (GenBank: AFRR00000000) [33]; the porcine lung isolates HN06, (GenBank: NC_017027.1) [26] and 3480, (GenBank: NC_017764.1); the bovine respiratory disease isolate 36950, (Genbank: NC_016808.1) [27] and the caprine isolate Anand1G, (GenBank: AFRS00000000). The uncorrected p method, which computes the proportion of positions at which two sequences differ, was used to build the distance matrix and determine distances between different strains. Phylogenetic trees were constructed with SplitsTree v4.11.3 [34] using the neighbour-joining method and uncorrected p distance matrices.
Results and Discussion
Genome sequencing of 12 P. multocida HS strains
Ten P. multocida strains, isolated from HS cases in buffalo (8 isolates) and cattle (2 isolates) from Pakistan and Thailand, and two Pakistani HS vaccine strains were used for this study (Table 1). Except for the Faisalabad isolate (Table 1), all of the isolates had previously been genotyped using MLST as ST122 [14]. Therefore, in order to identify if there were differences between the strains, we determined whole genome draft sequences of each strain. Genomic DNA was isolated from each strain and sequenced on an Illumina HiSeq 2000. Sequences reads were de novo assembled, resulting in between 32 (strains THF and V1) and 77 (strain Karachi) contigs of > 200 bp. The genomic features of the 12 sequenced strains and accession numbers are shown in Table 3. The predicted genome sizes ranged from 2.34 to 2.46 Mbp and the number of coding sequences (CDS) ranged between 2,082 (strains THA, THD and THF) and 2,179 (strain TX1). The GC content was highly conserved across all strains (40.31 to 40.41%).
Variant calls and genetics of capsule and LPS biosynthesis in the HS-causing strains
Mapping of the sequencing reads of the 12 Asian genomes to the complete genome of Pm70 identified between 16,443 and 16,513 CDS SNPs in the 12 genomes (Table 4). Furthermore, between 12 and 19 indels were also identified in the CDS of the 12 genomes (Table 4).
All HS-associated strains contained the type B capsular biosynthetic locus. Only three nucleotide changes were observed in the cap locus of all of the Asian strains compared to the cap locus of M1404 [35]. Two of the mutations were silent, while the third encoded a missense mutation (D50Y) within the putative glycosyltransferase EpsJ. An additional nucleotide change in a non-coding region was observed in the Pakistani strains (relative position: 1185274 in BUKK). The FASTA file of the cap locus of strain BUKK is provided as S1 File.
All HS-associated strains also contained the LPS outer core biosynthesis locus that is shared by strains belonging to Heddleston serovars 2 and 5 [9]. Only a single nucleotide change was observed across the LPS outer core biosynthesis loci of the 12 Asian strains when compared with the M1404 locus [9]. The FASTA file of the LPS biosynthesis locus of strain BUKK is provided as S2 File. The Heddleston 2 and 5 type strains can be differentiated serologically by the presence (serovar 5) or absence (serovar 2) of PEtn on the 2nd inner core heptose of the LPS [9]. Addition of this PEtn residue to the LPS is dependent on the presence of an intact lpt-3 gene (annotated as dcaA in Pm70) [9]. All HS strains contained a disrupted lpt-3 gene, with a nonsense mutation (relative position: 499751 in BUKK strain) identical to the mutation previously reported in M1404 [9]. Therefore, all of the HS strains analysed in this study are predicted to belong to LPS serovar 2. Indeed, all of the Pakistani and Thai strains had previously been classified using Heddleston serology as LPS serovar 2, except for THD which had been reported as LPS serovar 2,5 (P. Pathanasophon, personal communication, June, 2012). However, our genetic analyses would indicate that THD also belongs to LPS serovar 2, highlighting the advantage of molecular over serological typing methods.
It has been shown previously that most P. multocida strains produce two inner core LPS glycoforms designated A and B [30]. Production of these two glycoforms is dependent on the presence of two active heptosyltransferases, HptA and HptB. HptA is specific for inner core glycoform A and is responsible for the addition of the first heptose to a single phosphorylated 3-deoxy-D-manno-octulosonic acid (Kdo). For glycoform B assembly, unphosphorylated Kdo residues have a second Kdo residue added followed by the addition of the first heptose by a different heptosyltransferase, HptB [30]. The hptB gene in all strains included in this study was intact, indicating that HptB would be fully functional. The hptA gene was intact in all strains included in this study except for strain 3480 where it contained a premature stop codon (mutation at 2232375). Therefore, we would predict that strain 3480 is unable to add the first heptose to the glycoform A inner core. Interestingly, a hptA mutant constructed in the P. multocida fowl cholera isolate VP161 was significantly attenuated for virulence, predicted to be due to the presence of truncated glycoform A LPS on the surface of the cell [30]. However, in a further study it was found that growth of the hptA mutants in vivo selected for virulent strains with nonsense suppressor mutations in the Kdo phosphokinase gene, kdtA. This mutation prevented the phosphorylation of any Kdo residues, allowing all Kdo residues to be available for glycoform B LPS assembly [36].
Phylogeny of the tested strains
The relatedness of the different strains was determined by comparing all nucleotide changes at positions that were conserved across all of the comparison genomes (core genome single nucleotide polymorphisms; CG-SNPs). Firstly, all the P. multocida strains were compared together with the closely related species P. dagmatis, P. bettyae, Gallibacterium anatis and Mannheimia haemolytica. This analysis, using 789 CG-SNPs, clearly showed that all the P. multocida strains form a monophyletic group most closely related to P. dagmatis (Fig 1A). Secondly, comparison of only P. multocida strains using 7,892 CG-SNPs (Fig 1B), indicated that the P. multocida HS strains were very closely related and clearly separated from all of the other P. multocida strains. This finding is in contrast to previous analyses using a smaller number of strains that showed little or no correlation between phylogeny and serovar, disease type or host predilection [13]. However, there was still no clear correlation between strain relatedness and disease type other than for the HS strains. Indeed, the five fowl cholera isolates (Pm70, X73, VP161, P1059 and Anand1P) did not cluster separately from strains associated with other disease types. Finally, we compared just the HS strains sequenced in this study and M1404 using 722 CG-SNPs (Fig 1C). There was a clear separation between strain M1404 (the North American isolate) and the Thai and Pakistani strains, which were also clearly separated from each other. Therefore, the higher resolution provided by whole genome sequencing revealed a clear genetic relationship with geographic source which was not possible with MLST.
A. Relationship between Gallibacterium anatis, Mannheimia haemolytica, Pasteurella dagmatis, Pasteurella bettyae and the P. multocida strains Pm70, 36950, HN06, P3480, X73, VP161, Anand1C, Anand1B, Anand1P, Anand1G, P1059, P52VAC, VTCCBAA264, M1404 and the twelve Pakistani and Thai isolates. B. Relationship between the P. multocida strains. C. Relationship between the HS-associated P. multocida B:2 strain M1404 and the twelve Pakistani and Thai isolates. Phylogenetic relatedness for all comparisons was determined by analysis of only the single nucleotide polymorphisms found at conserved positions in all genomes of the comparison set (CG-SNPs); 789 shared positions for the tree in panel A, 7,829 shared positions for the tree in panel B and 722 shared positions for the tree in panel C. Trees were rendered with SplitsTree v4.11.3 [34]. The line segments above the trees with the number '0.01' indicate the branch length representing a genetic change of 0.01.
Core and pan genome predictions
We analysed the gene content of each of the 12 sequenced HS strains and compared these predictions with the coding sequences predicted for the four complete and annotated P. multocida genomes (36950, Pm70, 3480, HN06) and the M1404 draft genome. These analyses identified a shared set of 1,824 genes (core) and a pan genome of more than 2,700 genes (Fig 2). Furthermore, with the exception of the Thai isolates THA, THD and THF and the Pakistani isolates Pesh, Islm, V1, Faisal and ATTK, all other strains contained genes not found in any of the other strains used in the analysis (Fig 2). The unique genes in each of the four strains (BUKK, PVAcc, Karachi and Tx1) are provided in S1 Table. The genome from strains 36950 and 3480 contained 84 and 85 unique genes respectively; Pm70 contained 59 unique genes; strain TX1, 32 genes; strain Karachi, 11 genes and one unique gene was found in each genome of strains M1404, BUKK and PVAcc.
The central circle shows the number of genes common to all strains while the petals show the number of genes in addition to the core set, as well as the number of genes unique to each strain (in brackets). Abbreviated strain names are given outside each petal, strain details are given in Table 1. The B:2 HS-related strains are shaded as follows; M1404 is orange, Thai strains are blue and the Pakistani strains are green. Non-B:2 strains are shaded in grey.
Functional comparison of the core genes and strain-specific genes showed that the core genes are mainly responsible for inorganic ion transport and metabolism, energy production and metabolism, cell membrane biogenesis, ribosomal biogenesis, amino acid transport and metabolism, coenzyme transport and metabolism, carbohydrate transport and metabolism, signal transduction, transcription, translation, replication and repair, lipid metabolism, membrane transport and hypothetical proteins. The two strain-specific genes for PVAcc and BUKK encode hypothetical proteins. The eleven unique genes for the Karachi strain encode eight putative transposases and three hypothetical proteins. For strain TX1, the 32 unique genes encode putative conjugal transfer proteins, a DNA topoisomerase III, transcriptional regulator proteins, proteins involved in DNA replication, a cobalt ABC transport system, helicase and an endonuclease.
These comparative analyses also showed that all HS-associated strains included in the pan-genome analysis (M1404, Pakistani and Thai strains) share two large regions of unique sequence compared to the other complete genomes. The first region is approximately 34 kb in length (region 3 in Fig 3) while the second region is approximately 15 kb in length (region 4 in Fig 3) (Table 5). Furthermore, there are several dispersed genes uniquely present in all of the HS strains. Overall the HS strains share 96 genes that are absent from the other genomes analysed, including the capsule biosynthesis locus, present in all strains belonging to capsular serogroup B (Fig 3). These genes unique to the HS strains are provided in S1 Table. In addition, the twelve Asian HS strains share an approximately 44 kb region (region 2 in Fig 3) (Table 5) that is absent from the American HS strain M1404; this region contains 59 genes that encode mostly proteins with no significant similarity to proteins of known function. The seven Pakistani strains (ATTK, Faisal, Islm, Karachi, Pesh, PVAcc and V1) share an approximately 50 kb region containing 39 genes that is not present in the other genomes (region 1 in Fig 3) (Table 5). Genes encoded in region 1 encode mostly phage elements as well as hypothetical proteins. Additionally, strains 36950, TX1 and BUKK share 35 unique genes, encoding elements with similarity to the integrative conjugative element (ICEPmu1) of 36950 (Fig 4). There are also 42 genes shared by TX1 and BUKK. TX1 is the only Asian HS strain with a large number (32) of unique genes and these predominantly encode phage elements and hypothetical proteins. The three Thai HS isolates have just a single unique gene encoding an abortive infection-(Abi-) like protein of 226 amino acids. Abi-like genes are found in various bacterial species, and encode proteins involved in bacteriophage resistance [37,38].
The three inner rings show the DNA size, GC content and GC skew of the reference genome (PVAcc strain). The 17 outer rings show regions of the comparison genomes that match the reference genome PVAcc and in the order (inside to outside) 3480, 36950, HN06, Pm70, M1404, BUKK, Tx1, THA, THD, THF, ATTK, Faisal, Islm, Pesh, PVAcc, V1 and Karachi. Regions 1–4 on the outside identify particular regions of difference between the strains that are potential prophage elements. The position of the type B capsule locus is also noted between regions 2 and 3. This figure was drawn using BLAST ring image generator (BRIG) [39].
The three inner rings show the DNA size, GC content and GC skew of the reference element (ICEPmu1). The five outer rings show regions of the comparison genomes that match the reference ICEPmu1. Figure was drawn using BRIG [39].
Phage identification
All of the genomes from known HS-associated strains were analysed for the presence of phage elements using PHAST [28]. This analysis identified four regions corresponding to putative temperate phage elements. The genomic locations of these four regions in PVAcc strain, as an example, are presented in Table 5. These regions correspond with the main genetic differences identified between different groups of strains (regions 1–4; Fig 3). Regions 3 and 4 were identified as intact phages and have been reported previously [13]. Region 2 was identified as an incomplete phage and region 1 as a questionable phage (Fig 3). The questionable phage identified in region 1 is present in the seven Pakistani strains, the incomplete phage identified in region 2 is shared by all the Asian strains (relative position BUKK_04735–04880) and the intact phages identified in region 3 and 4 are shared by all the HS strains (relative position BUKK_07250–07540). The phage elements identified in regions 1, 2, 3 and 4 are situated at tRNALeu, tRNAMet, tRNASer and tRNAMet genes respectively. This correlates with the previous reports that the F108 phage and the lysogenic phage carrying the PMT toxin, integrate into the t33tRNALeu and the t3tRNALeu genes respectively [13,40]. Temperate phages may contain important virulence genes [41]. Indeed, as noted above the P. multocida PMT toxin is the primary virulence factor for porcine atrophic rhinitis and is carried on a lysogenic bacteriophage [42]. Therefore, the presence of different phage elements in different sets of strains may impact on the virulence of these strains. Further studies should assess the impact of these different gene sets on virulence.
A P. multocida HS-specific diagnostic PCR has been developed previously [15]. We searched for the DNA sequence recognized by this PCR in all genomes and identified it within the putative intact prophage within region 3 (Fig 3), Thus, this sequence is indeed present in all of the HS strains analysed in this study and is not present in any of the non HS-associated genomes analysed.
Identification of antibiotic resistance genes and characterization of ICEPmu2
The 13 HS-associated genomes were analysed for the presence of acquired antimicrobial resistance genes using the ResFinder tool [29]. It should be noted that ResFinder searches only for acquired resistance genes and not for mutations in chromosomally-encoded genes that lead to antibiotic resistance. Acquired antimicrobial resistance genes were identified only in the BUKK and Tx1 Pakistani isolates. These included three aminoglycoside resistance genes (strA, strB and aph(3’)-lc), one beta lactamase gene (blaTEM-1B), one chloramphenicol resistance gene (catA2), one sulphonamide resistance gene (sul2) and one tetracycline resistance gene (tet(H)). Thus, these two strains should be resistant to streptomycin, kanamycin/neomycin, beta-lactams, chloramphenicol, sulphonamides and tetracycline. Indeed, this correlates with the clinical data on these isolates as the use of beta lactam antibiotics for infections involving BUKK and TX1 strains has been avoided due to the beta lactamase activity exhibited by these strains (E. Nawaz, personal communication, June, 2012).
All of the identified antibiotic resistance genes were clustered in two regions (relative position BUKK_06700–06910 and 11175–11375). However, it is likely that these regions are colinear within each genome but, due to contig breaks, have been separated in each draft genome assembly. We propose that this region encodes an ICE as many of the genes in this region encode proteins with shared identity to those in the ICEPmu1 of strain 36950 [27]. Using 50% identity and 50% gene coverage as the minimum criteria for a match, a total of 77 genes were identified in this region, all of which were shared by the BUKK and TX1 HS strains. Of these, 35 encoded proteins that shared a significant level of identity with proteins in the ICEPmu1 region of the multi-drug resistant strain 36950 and a further 16 encoded proteins with lower levels of identity (30 to 50%) to proteins in ICEPmu1 (Fig 4). Thus, it is likely that the BUKK and TX1 strains both contain an ICE. The genes with shared identity to the ICEPmu1 genes are located in a number of colinear groups within this ICE. These groups are interspersed with the genes identified as unique to strains BUKK and TX1 (Fig 4). The end of one of these colinear groups is flanked by tRNALeu. The equivalent region in strain 36950 represents the right end of the ICEPmu1 element (Pmu_03540–03610) and encodes a number of proteins involved in DNA replication, including a single-stranded DNA-binding protein, an ATPase involved in chromosome partitioning, a DnaB-like helicase and a ParB family protein with a predicted DNA nuclease domain. This set of genes has been reported as the most conserved region among diverse proteobacterial ICE [27,43].
Within the putative ICE identified in BUKK and TX1, 47 genes encoded proteins with predicted functions. These included proteins predicted to be involved in ICE mobility, including excision/integration and conjugative transfer. A putative phage integrase was identified (BUKK_06905) with similarity to two integrases found in the ICEPmu1 within strain 36950 (identities of 51% and 57% with the first and second integrases, respectively). BUKK_06905 shared significant similarity with tyrosine recombinases of the Xer family, which mediate integration via site-specific recombination. A gene encoding a putative relaxase (BUKK_06900) was identified downstream of this integrase gene; a similar organisation is found in the ICEPmu1 [27]. Proteins necessary for conjugative transfer were also present, including proteins predicted to be necessary for the formation of a type IV pilus (BUKK_11185 which shows 53% similarity with Pmu_03230), TraD (BUKK_11205 which shows 73% similarity with Pmu_03190), TraG (BUKK_11275 which shows 60% similarity with Pmu_03040) and TraC (BUKK_11210 which shows 67% similarity with Pmu_03070) [27]. Moreover, a gene encoding a protein with a lysozyme-like domain (BUKK_11195 which shows 52% similarity with Pmu_03210) and a putative DNA topoisomerase III (BUKK_06775) (66% similarity with Pmu_03290 in strain 36950) were also identified.
While the resistance genes strA (BUKK_06815) and strB (BUKK_06820), aph(3’)-lc (BUKK_11370), sul2 (BUKK_06810) and tetR-tet(H) (BUKK_11315 and BUKK_11320) are present in strains BUKK, TX1 and 36950 [27,44], the blaTEM-1B (BUKK_06875) and catA2 (BUKK_11355) genes are unique to strains BUKK and TX1. However, strain 36950 contains more resistance genes than strains BUKK and TX1 as it also contains resistance genes for streptomycin/spectinomycin (aadA25), gentamicin (aadB), kanamycin/neomycin (aphA1), chloramphenicol/florfenicol (floR), tilmicosin/clindamycin (erm(42)) and tilmicosin/tulathromycin (msr(E)-mph(E)) [44]. These differences indicate that the putative ICE present in strains BUKK and TX1 is not identical to ICEPmu1 in strain 36950 and should therefore be designated ICEPmu2 as the second ICE discovered in P. multocida. Further work should aim to close the contig breaks in this region to confirm this is a single element. In addition it would be of great interest to investigate the mobility of ICEPmu2.
Conclusion
In conclusion, we have shown that HS-associated P. multocida strains belonging to capsular serogroup B form a very closely related group, but are distinguishable using whole genome analyses. We identified 96 genes unique to the HS-associated strains and future characterization of these genes should elucidate the roles they play in disease pathogenesis, virulence and host specificity. Selected genes from this group will be excellent candidates for the development of a rapid diagnostic test for HS. The putative integrative conjugative element (ICE) identified in two Pakistani isolates should be further analysed to determine its mobility and relatedness to ICEPmu1 of strain 36950.
Supporting Information
S1 File. Nucleotide sequence of the strain BUKK capsule biosynthesis locus in FASTA format.
https://doi.org/10.1371/journal.pone.0130296.s001
(TXT)
S2 File. Nucleotide sequence of the strain BUKK LPS biosynthesis locus in FASTA format.
https://doi.org/10.1371/journal.pone.0130296.s002
(TXT)
S1 Table. Genes unique to the HS strains, TX1 strain, Karachi strain, BUKK strain and PVacc strain.
https://doi.org/10.1371/journal.pone.0130296.s003
(XLSX)
Acknowledgments
A special gratitude needs to go to Prof. John Edwards, Prof. Stan Fenwick, Dr Tim Hyndman and Dr Kirsty Townsend who helped at various stages of this project. To Prof. M. Iqbal Choudhary, the head of the International Center for Chemical and Biological Sciences (ICCBS), University of Karachi, Pakistan, and to Dr Kamran Azim and Dr Muhammed A. Mesaik, thanks for facilitating the work done in Pakistan. Thanks to the National Veterinary Laboratories (NVL) in Islamabad, Pakistan, for providing the Pakistani P. multocida isolates and the National Institute of Animal Health (NIAH) at the Department of Livestock Development in Bangkok, Thailand, for providing the Thai P. multocida isolates. Thanks also to the Veterinary Research Institute (VRI) in Lahore, Pakistan, for providing a vaccine strain and the Beijing Genomics Institute (BGI), China, for whole genome sequencing.
Author Contributions
Conceived and designed the experiments: AM JB MB. Performed the experiments: AM TS SG. Analyzed the data: AM TS SG BA MH JB MB. Wrote the paper: AM TS SG BA MH JB MB.
References
- 1.
Rimler RB, Rhoades KR (1989) Pasteurella multocida. In: Adlam C, Rutter JM, editors. Pasteurella and Pasteurellosis: Academic press, London. pp. 37–73.
- 2.
De Alwis MCL (1999) Haemorrhagic septicaemia. Canberra: Australian Centre for International Agricultural Research. pp. 11–24.
- 3. Carter GR. A haemagglutination test for identification of serological types. Am J Vet Res. 1955; 16: 481–484. pmid:13238757
- 4. Townsend KM, Boyce JD, Chung JY, Frost AJ, Adler B. Genetic organization of Pasteurella multocida cap loci and development of a multiplex capsular PCR typing system. J Clin Microbiol. 2001; 39: 924–929. pmid:11230405
- 5. Heddleston KL, Gallagher JE, Rebers PA. Fowl cholera: Gel diffusion precipitin test for serotyping Pasteruella multocida from avian species. Avian Dis. 1972; 16: 925–936. pmid:4628021
- 6. Benkirane A, De Alwis MCL. Haemorrhagic septicaemia, its significance, prevention and control in Asia. Vet Med (Praha). 2002; 47: 234–240.
- 7.
Neramitmansook P (1993) Thailand. In: Patten BE, Spencer TL, Johnson RB, Hoffmann D, Lehane L, editors. Pasteurellosis in production animals, an international workshop sponsored by ACIAR held at Bali, Indonesia, 10–13 August 1992. Canberra: ACIAR Proceedings No. 43. pp. 234–237.
- 8. Dziva F, Muhairwa AP, Christensen H, Bisgaard M. Diagnostic and typing options for investigating diseases associated with Pasteurella multocida. Vet Microbiol. 2008; 128: 1–22. pmid:18061377
- 9. St Michael F, Harper M, Parnas H, John M, Stupak J,Vinogradov E, et al. Structural and genetic basis for the serological differentiation of Pasteurella multocida Heddleston serotypes 2 and 5. J Bacteriol. 2009; 191: 6950–6959. pmid:19767423
- 10. May BJ, Zhang Q, Li LL, Paustian ML, Whittam TS, Kapur V. Complete genomic sequence of Pasteurella multocida,Pm70. PNAS. 2001; 98: 3460–3465. pmid:11248100
- 11.
NCBI Genome website. Available: http://www.ncbi.nlm.nih.gov/genome/genomes/912. Accessed 17 June 2015.
- 12. Wilson BA, Ho M. Pasteurella multocida: From zoonosis to cellular microbiology. Clin Microbiol Rev. 2013; 26: 631–655. pmid:23824375
- 13.
Boyce JD, Seemann T, Adler B, Harper M (2012) Pathogenomics of Pasteurella multocida. In: Aktories K, Orth JHC, Adler B, editors. Pasteurella multocida: Molecular Biology, Toxins and Infection: Springer Berlin Heidelberg. pp. 23–38.
- 14. Moustafa AM, Bennett MD, Edwards J, Azim K, Mesaik MA, Choudhary MI, et al. Molecular typing of haemorrhagic septicaemia-associated Pasteurella multocida isolates from Pakistan and Thailand using multilocus sequence typing and pulsed-field gel electrophoresis. Res Vet Sci. 2013; 95: 986–990. pmid:23916592
- 15. Townsend KM, Frost AJ, Lee CW, Papadimitriou JM, Dawkins HJS. Development of PCR assays for species- and type-specific identification of Pasteurella multocida isolates. J Clin Microbiol. 1998; 36: 1096–1100. pmid:9542944
- 16. Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012; 19: 455–477. pmid:22506599
- 17. Zerbino DR, Birney E. Velvet: Algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 2008; 18: 821–829. pmid:18349386
- 18. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990; 215: 403–410. pmid:2231712
- 19.
NCBI Blast website. Available: http://blast.ncbi.nlm.nih.gov/Blast.cgi. Accessed 17 June 2015.
- 20. Gurevich A, Saveliev V, Vyahhi N, Tesler G. QUAST: quality assessment tool for genome assemblies. Bioinformatics. 2013; 29: 1072–1075. pmid:23422339
- 21. Rissman AI, Mau B, Biehl BS, Darling AE, Glasner JD,Perna NT, et al. Reordering contigs of draft genomes using the Mauve Aligner. Bioinformatics. 2009; 25: 2071–2073. pmid:19515959
- 22. Angiuoli SV, Gussman A, Klimke W, Cochrane G, Field D, Garrity G, et al. Toward an online repository of standard operating procedures (SOPs) for (meta)genomic annotation. OMICS. 2008; 12: 137–141. pmid:18416670
- 23.
NCBI Genbank website. Available: http://www.ncbi.nlm.nih.gov/genbank. Accessed 17 June 2015.
- 24. Benson DA, Cavanaugh M, Clark K, Karsch-Mizrachi I, Lipman DJ, Ostell J, et al. GenBank. Nucleic Acids Res. 2013; 41: D36–D42. pmid:23193287
- 25. Darling AE, Mau B, Perna NT. progressiveMauve: Multiple genome alignment with gene gain, loss and rearrangement. PLoS One. 2010; 5: e11147. pmid:20593022
- 26. Liu W, Yang M, Xu Z, Zheng H, Liang W,Zhou R, et al. Complete genome sequence of Pasteurella multocida HN06, a toxigenic strain of serogroup D. J Bacteriol. 2012; 194: 3292–3293. pmid:22628516
- 27. Michael GB, Kadlec K, Sweeney MT, Brzuszkiewicz E, Liesegang H,Daniel R, et al. ICEPmu1, an integrative conjugative element (ICE) of Pasteurella multocida: Structure and transfer. J Antimicrob Chemother. 2012; 67: 91–100. pmid:22001176
- 28. Zhou Y, Liang Y, Lynch KH, Dennis JJ, Wishart DS. PHAST: A fast PHAge Search Tool. Nucleic Acids Res. 2011.
- 29. Zankari E, Hasman H, Cosentino S, Vestergaard M, Rasmussen S, Lund O, et al. Identification of acquired antimicrobial resistance genes. J Antimicrob Chemother. 2012; 67: 2640–2644. pmid:22782487
- 30. Harper M, Boyce JD, Cox AD, St. Michael F, Wilkie IW, Blackall PJ, et al. Pasteurella multocida expresses two lipopolysaccharide glycoforms simultaneously, but only a single form is required for virulence: Identification of two acceptor-specific heptosyl I transferases. Infect Immun. 2007; 75: 3885–3893. pmid:17517879
- 31. Vaid RK, Shanmugasundaram K, Boora A, Bera BC, Shukla BN, Anand T, et al. Draft genome sequence of Pasteurella multocida subsp. multocida B:2 strain VTCCBAA264 isolated from Bubalus bubalis in North India. Genome Announcements. 2014; 2.
- 32. Abrahante JE, Johnson TJ, Hunter SS, Maheswaran SK, Hauglund MJ, Bayles DO, et al. Draft genome sequences of two virulent serotypes of avian Pasteurella multocida. Genome Announcements. 2013; 1.
- 33. Ahir VB, Roy A, Jhala MK, Bhanderi BB, Mathakiya RA, Batt VD, et al. Genome sequence of Pasteurella multocida subsp. gallicida Anand1_poultry. J Bacteriol. 2011; 193: 5604. pmid:21914901
- 34. Huson DH, Bryant D. Application of Phylogenetic Networks in Evolutionary Studies. Mol Biol Evol. 2006; 23: 254–267. pmid:16221896
- 35. Boyce J, Chung JY, Adler B. Genetic organisation of the capsule biosynthetic locus of Pasteurella multocida M1404 (B:2). Vet Microbiol. 2000; 72: 121–134. pmid:10699509
- 36. Harper M, Cox AD, St Michael F, Ford M, Wilkie I, Adler B, et al. Natural selection in the chicken host identifies 3-deoxy-D-manno-octulosonic acid kinase residues essential for phosphorylation of Pasteurella multocida lipopolysaccharide. Infect Immun. 2010; 78: 3669–3677. pmid:20566690
- 37. Anba J, Bidnenko E, Hillier A, Ehrlich D, Chopin MC. Characterization of the lactococcal abiD1 gene coding for phage abortive infection. J Bacteriol. 1995; 177: 3818–3823. pmid:7601848
- 38. Garvey P, Fitzgerald GF, Hill C. Cloning and DNA sequence analysis of two abortive infection phage resistance determinants from the lactococcal plasmid pNP40. Appl Environ Microbiol. 1995; 61: 4321–4328. pmid:8534099
- 39. Alikhan N-F, Petty N, Ben Zakour N, Beatson S. BLAST Ring Image Generator (BRIG): Simple prokaryote genome comparisons. BMC Genomics. 2011; 12: 402. pmid:21824423
- 40. Campoy S, Aranda J, Àlvarez G, Barbé J, Llagostera M. Isolation and sequencing of a temperate transducing phage for Pasteurella multocida. Appl Environ Microbiol. 2006; 72: 3154–3160. pmid:16672452
- 41. Boyd EF, Brüssow H. Common themes among bacteriophage-encoded virulence factors and diversity among the bacteriophages involved. Trends Microbiol. 2002; 10: 521–529. pmid:12419617
- 42. Pullinger GD, Bevir T, Lax AJ. The Pasteurella multocida toxin is encoded within a lysogenic bacteriophage. Mol Microbiol. 2004; 51: 255–269. pmid:14651626
- 43. Mohd-Zain Z, Turner SL, Cerdeño-Tárraga AM, Lilley AK, Inzana TJ,Duncan AJ, et al. Transferable antibiotic resistance elements in Haemophilus influenzae share a common evolutionary origin with a diverse family of syntenic genomic islands. J Bacteriol. 2004; 186: 8114–8122. pmid:15547285
- 44. Michael GB, Kadlec K, Sweeney MT, Brzuszkiewicz E, Liesegang H, Daniel R, et al. ICEPmu1, an integrative conjugative element (ICE) of Pasteurella multocida: Analysis of the regions that comprise 12 antimicrobial resistance genes. J Antimicrob Chemother. 2012; 67: 84–90. pmid:22001175