Genome sequence analysis of the Indian strain Mannheimia haemolytica serotype A2 from ovine pneumonic pasteurellosis

Abstract

Mannheimia haemolytica is a leading causative agent of pasteurellosis in ruminants. Genome of M. haemolytica strains from different hosts has been sequenced worldwide to understand its pathogenesis. There are only few reports on the isolation of M. haemolytica in India with limited information on its molecular characteristics. The present study focuses on genome sequence analysis of a M. haemolytica strain isolated from pneumonic sheep. Mannheimia haemolytica A2 strain NIVEDI/MH/1 was isolated and identified by species and serotype-specific PCRs. Whole genome sequencing was performed using the Ion Torrent Personal Genome Machine. A comparative genomic analysis was performed to understand the virulence determinants of the Indian strain and its phylogenetic relationship with other global strains. Sequence data revealed a draft genome of 2,211,426 bp size with 41.3% GC content, assembled into 17 contigs, and contained 2379 genes. Five genomic islands identified in the genome showed high sequence identity with other respiratory pathogens of the Pasteurellaceae family. Phylogenetic analysis showed M. haemolytica A2 NIVEDI/MH/1 is very close to a M. haemolytica A2 strain from pneumonic calf. Further, the analysis revealed the presence of virulence, metal-, and multidrug resistance genes needed for pathogenesis and survival of the bacteria during infection. Also, we identified the presence of type I-C and type II-C of CRISPR-Cas arrays in the present sequenced genome. The study emphasizes the role of M. haemolytica in respiratory infections of ruminants in the Indian subcontinent and indicates the role of vertical and horizontal gene pools in pathogenicity and survivability of the bacteria.

Introduction

Mannheimia haemolytica is a Gram-negative, weakly hemolytic coccobacillus (Highlander 2001). It is an important etiological agent of pasteurellosis, also known as respiratory mannheimiosis (Mohamed and Abdelsalam 2008) and mastitis in ruminants (Omaleki et al. 2011, 2016). Based on the capsular antigens, M. haemolytica is classified into 12 serotypes (A1, A2, A5-A9, A12- A14, A16, and A17) (Klima et al. 2017). Of these, the serotypes A1 and A6, responsible for the bovine respiratory disease (BRD) in cattle, are extensively studied, whereas the serotypes A2 and A7, causing pneumonic pasteurellosis in small ruminants, are least characterized (Highlander 2001). Stress due to unfavorable environmental conditions, animal transportation, and prior viral and bacterial infections act as predisposing factors in triggering this opportunistic pathogen to cause disease outbreaks (Timsit et al. 2013). The molecular mechanisms involved in the conversion of this bacteria from commensal to pathogen are still unclear (Gioia et al. 2006).

The virulence factors determine the pathogenicity of an organism (Mohamed and Abdelsalam 2008). These virulence factors of M. haemolytica include ruminant-specific leukotoxin, anti-phagocytic capsule, lipopolysaccharide, iron-regulated outer membrane proteins, lipoproteins, sialoglycoprotease, neuraminidase, and two immunoglobulin proteases (Highlander 2001). In addition, the emergence of multi-antibiotic resistance to beta-lactams, tetracycline, streptomycin, sulfonamides, macrolides, and sulfamethazine in M. haemolytica has been reported (Hörmansdorfer and Bauer 1996).

Genomes of different serotypes of M. haemolytica were sequenced to identify novel targets for diagnostics and vaccine development (Klima et al. 2016), and to understand the phylogenetic relationship with other strains reported across the globe (Moustafa et al. 2015). Genome sequence analysis of M. haemolytica A1, strain ATCC BAA-410, from a lung of a calf that died of BRD identified the presence of virulence, natural competence, accessory, and transcription regulation genes and revealed widespread genes conserved in M. haemolytica as with other members of the Pasteurellaceae family. The study also explained the similarity of M. haemolytica virulence proteins with Neisseria meningitidis, suggesting a convergent evolution (Gioia et al. 2006). Lawrence et al. (2010b) reported the presence of phage regions and pseudogenes in the genome of M. haemolytica serotype A2, indicating pathogenic islands as a consequence of the horizontal gene transfer (HGT)). HGT is an important factor in the emergence of the pathogen by the evolution of the novel virulence genes and the spread of antimicrobial and metal resistance genes (Keen 2012). A previous study has substantiated the role of Mannheimia genomic islands (GIs) and HGT in the expression and regulation of virulence factors (Rao and Jayakumar 2017). Antibiotic and metal resistance genes were found to be associated with plasmids, chromosome (Watts et al. 1994), and integrative conjugative elements (ICEs) (Eidam et al. 2015). M. haemolytica genome also contains the CRISPR-Cas systems (clustered regularly interspaced short palindromic repeats/CRISPR-associated proteins) (Klima et al. 2016), which act as an adaptive immune system for bacterial protection against invading genetic elements like phages and plasmids (Makarova et al. 2011).

There are fewer reports on the prevalence of M. haemolytica in India (Sharma et al. 2011; Dar et al. 2012; Kumar et al. 2015). However, no reports are available on the genome sequencing and the molecular characterization of M. haemolytica isolated from India. Hence, the present study was focused on the whole genome sequence analysis of M. haemolytica isolated from a pneumonic sheep to unravel the genes involved in the virulence and its phylogenetic relatedness to the global isolates.

Materials and methods

Sheep herd

The sheep samples were collected from a private herd at Chikkajala village, Urban Bengaluru, Karnataka. The herd was maintained with a semi-intensive rearing method. The animals were not vaccinated and a sudden death of six sheep due to the respiratory infection was reported. The sample was collected from the sheep having the typical case of uncomplicated pneumonia with symptoms including dyspnea, pyrexia, and mucopurulent nasal discharge.

Isolation of M. haemolytica and extraction of genomic DNA

M. haemolytica was isolated from the nasopharynx of the sheep (Ovis aries) with a history of pneumonia. The isolate was tentatively identified as M. haemolytica based on conventional culture characteristics including colony morphology and standard biochemical tests (Marru et al. 2013), and confirmed by species and capsule-specific PCR assays (Alexander et al. 2008; Klima et al. 2017). The M. haemolytica A2 isolate was designated as strain NIVEDI/MH/1. The strain was initially cultured on tryptic soy agar (TSA) supplemented with 5–7% sheep blood at 37 °C for 24 h and subsequently subcultured on brain heart infusion agar (BHI) for genomic DNA extraction. DNA was extracted using the DNeasy kit as per the manufacturer’s protocol (Qiagen, Hilden, Germany). The quantity and quality of the extracted DNA were determined by using NanoDrop2000 (Thermo Scientific, Waltham, USA) and by agarose gel (0.8%) electrophoresis, respectively.

Genome sequencing, assembly, and annotation

The whole genome of M. haemolytica strain A2 NIVEDI/MH/1 was sequenced using the Ion Torrent Personal Genome Machine (Life Technologies, Carlsbad, CA). The de novo assembly was performed using MIRA (Mimicking Intelligent Read Assembly) version 3.9.18 (Chevreux et al. 1999). For the genome annotation, the assembled sequence of M. haemolytica A2 NIVEDI/MH/1 was submitted to the NCBI Prokaryotic Genomes Annotation Pipeline (NCBI PGAP) (http://www.ncbi.nlm.nih.gov/ genome/ annotation prok/ process/). Further, annotation was also performed by using Rapid Annotations Subsystems Technology (RAST) (Aziz et al. 2008). The rRNAs and tRNAs genes were predicted using RNAmmer (Lagesen et al. 2007) and tRNAscan-SE 1.21 (Lowe and Eddy 1997), respectively. Protein-coding sequences were analyzed by the Cluster of Orthologous Group (COG) database (Tatusov et al. 2000) on WebMGA (Wu et al. 2011). Graphical view of the genome was generated using Artemis and DNAplotter. The total number of contigs, their predicted coding sequences in forward and reverse direction, rRNA, and tRNA were plotted.

16S rRNA and genome-based phylogenetic analysis

The 16S rRNA gene from M. haemolytica A2 NIVEDI/MH/1 genome was identified using RNAmmer. The 16S rRNA genes of other genome sequences were retrieved from the NCBI RefSeq database (Supplementary Table S1). A phylogenetic tree based on the 16S rRNA genes was constructed using Molecular Evolutionary Genetics Analysis (MEGA) version 7 (Tamura et al. 2013) by the maximum likelihood method with the default bootstrap value as 500. The whole genome sequences of 62 M. haemolytica, 3 Mannheimia varigena, and one each of Mannheimia granulomatis, Mannheimia massilioguelmaensis, Mannheimia succiniciproducens, and Mannheimia sp. were retrieved from NCBI (Supplementary Table S1). These genome sequences were compared with the genome sequence of M. haemolytica A2 NIVEDI/MH/1. The Composition Vector (CV) method (Qi et al. 2004) was used to generate the distance matrix. The whole genome-based phylogenetic tree was constructed using the Neighbor-Joining (NJ) method and visualized using MEGA version 7.

Identification and prediction of virulent genes and genomic islands

The virulence factors were predicted using the Virulence Factor Database (VFDB) (Chen et al. 2016). The Genomic Islands (GIs) were envisaged using Island viewer 4, a computational tool integrated with four different genomic island prediction methods such as IslandPath-DIMOB, SIGI-HMM, Islander, and IslandPick. The IslandPath-DIMOB is based on several common characteristics of GIs such as abnormal sequence composition and presence of mobility genes (Hsiao et al. 2003). The SIGI-HMM is based on codon usage bias with the Hidden Markov Model approach (Waack et al. 2006). The Islander is based on mechanistic consequences of their typical site-specific integration into tRNA/tmRNA genes (Hudson et al. 2014), and the IslandPick is based on a comparative genomic approach (Langille et al. 2008). The Panseq server (Laing et al. 2010) was used to find the novel regions in the predicted GIs of M. haemolytica A2 NIVEDI/MH/1 by comparing with other existing genome sequences.

Prediction of the CRISPR-Cas system

The CRISPR-Cas system was predicted with the use of web tools CRISPRFinder (Grissa et al. 2007) and CRISPRone (Zhang and Ye 2017). CRISPRFinder, a web service tool, was used to detect CRISPR-Cas array in the sequence and define the direct repeats and spacers. CRISPRone was used for the annotation of the predicted CRISPR-Cas system to detect Cas genes and type of the CRISPR-Cas system present in the M. haemolytica A2 NIVEDI/MH/1 genome.

Results and discussion

Mannheimiosis is an emerging infectious disease of small ruminants in India. This study was focused on the isolation, identification, and whole genome sequencing of the M. haemolytica serotype A2 strain from a pneumonic sheep aimed to unravel its genetic background circulating in Indian terrain. The whole genome sequencing of M. haemolytica A2 NIVEDI/MH/1 using Ion Torrent PGM yielded a total of 1,556,783 reads with an average read length of 167 bp. The total sequence output was ~ 260 Mbp, which is approximately 189-fold coverage. The de novo assembly using MIRA yielded 17 contigs, and the largest contig was 586,935 bp long. The draft genome of M. haemolytica A2 NIVEDI/MH/1 was 2,211,426 bp with 41.32% GC content. Genome annotation by RAST predicted a total of 2379 protein-coding regions. Of these, 2023 genes were predicted with functions, and 356 genes were annotated as hypothetical proteins. RNAmmer and tRNAscan-SE 1.21 predicted a total of 17 rRNA and 58 tRNA genes, respectively (Fig. 1; Table 1). The COG category of the predicted genes is shown in Table 2.

Fig. 1

Circular representation of M. haemolytica A2 NIVEDI/MH/1 genome. From the outer circle to the inner circle, black indicates the total base pair, red indicates forward strand, blue indicates reverse strand, gray indicates a total number of contigs, pink indicates tRNA, green indicates rRNA, gray and light blue indicate GCskew, and green and pink indicate GCplot

Table 1 Genome statistics of M. haemolytica A2 NIVEDI/MH/1

Table 2 Number of genes associated with general COG functional categories

The 16S rRNA gene sequence can be used for the species identification and phylogeny analysis (Clarridge 2004). However, the whole genome sequences are generally preferred to study the phylogenetic relatedness on the basis of geographical origin, serotypes, and host predilection (Moustafa et al. 2015). The 16S rRNA gene-based phylogenetic tree confirmed that the strain NIVEDI/MH/1 belongs to the species M. haemolytica as it was clustered with other strains of M. haemolytica in the same clade (Fig. 2). The whole genome-based phylogenetic tree constructed using 70 Mannheimia strains showed that the strain M. haemolytica A2 NIVEDI/MH/1 was the closest neighbor of M. haemolytica D35 serotype A2 (accession no. AUNK00000000), which was isolated from a pneumonic calf lung (Hauglund et al. 2015). Similarly, the M. haemolytica serotype A2 ovine strain from a pneumonic lung (accession no. ACZX00000000) (Lawrence et al. 2010a) was closely related to the M. haemolytica A2 NIVEDI/MH/1 strain and clustered together in the same clade. The phylogenetic tree based on both the 16S rRNA gene and whole genome sequence showed the similarity of M. granulomatis with the M. haemolytica A1/A6 PLK10 strain (Supplementary Fig. S1) and three M. varigena strains forming a separate clade but clustering with Mannheimia sp.

Fig. 2

16S rRNA gene-based phylogenetic analysis of M. haemolytica A2 NIVEDI/MH/1 and other Mannheimia sp. genome retrieved from the NCBI RefSeq database (access date: 2/05/17) using the Maximum likelihood method in the MEGA software version 7. The analysis revealed that M. haemolytica A2 NIVEDI/MH/1 isolate (in red) clustered to the same clade with other M. haemolytica

A total of 33 virulent genes were predicted by the VFDB database (Supplementary File S2). Of these, lpxA (ODQ38859), lpxH (ODQ37174), lpxA (ODQ38611), lpxB (ODQ35950), rfaE (ODQ38829), wecA(ODQ36883), waaQ (ODQ37174), and lsgE (ODQ37931) are needed for the biosynthesis of lipopolysaccharide (LPS). Previous studies have reported that all these genes are conserved among the members of the Pasteurellaceae family (Xu et al. 2008; Peng et al. 2016). LPS is an endotoxin responsible for the stimulation of IL-1beta, IL-8, and neutrophil influx leading to inflammation and damage of bovine pulmonary endothelial cells (Paulsen et al. 1989; Yoo et al. 1995; Lafleur et al. 1998). The gene kdkA (ODQ37905) is needed for the LPS glycoform generation (Peng et al. 2016). According to a previous study, LPS enhances the cytolytic activity of leukotoxin, (repeats in-toxin (RTX toxin)) of M. haemolytica (Lafleur et al. 1998). RAST annotation identified 28 genes are having virulence potential (Supplementary File S3). NCBI PGAP and RAST annotations identified the leukotoxin-activating lysine-acyltransferase LktC (ODQ39212) needed for the posttranslational acylation of the leukotoxin protein encoded by LktA in the LktCABD operon (Highlander et al. 1989).

The presence of NanA (ODQ39248), NanK (ODQ39247), and a hypothetical protein (ODQ37835) having similarity with sialidase (neuraminidase) are needed for the sialometabolism that helps bacteria to colonize, persist, and cause disease in mammalian species (Sakarya and Oncü 2003; Steenbergen et al. 2005; Severi et al. 2007) were identified. The gene for UDP-N-acetylglucosamine 2-epimerase (ODQ37932.1), a rate-limiting enzyme in sialometabolism (Klima et al. 2016), was identified in M. haemolytica A2 NIVEDI/MH/1. Also, adhesin (ODQ39272) and a type IV pilus biogenesis/stability protein PilW (ODQ38480) required for adhesion and motility as noticed in Haemophilus influenzae, Pseudomonas aeruginosa, and Neisseria sp. were identified (Bakaletz et al. 2005). The genome contained the TonB transport system, TonB-ExbB-ExbD (ODQ38608, ODQ38777, and ODQ38778) (Peng et al. 2016), ATP-binding iron ABC transporter (ODQ38724, ODQ38726, and ODQ38478), heme utilization proteins (HutZ and Hut W) (ODQ38471, ODQ38473), and iron-binding protein (sitA) (ODQ38479) needed for the iron acquisition. It also contained a ferric uptake regulator (Fur family) (ODQ38822), a transcriptional repressor required for regulation of iron acquisition genes. Other putative virulence factors are superoxide dismutase (ODQ38739), virulence protein (ODQ38476), antitoxin, component of a toxin-antitoxin (TA) (ODQ38475), lipoprotein localization factor LolB (ODQ37719), S-ribosyl homocysteine lyase, Lux S (ODQ38823), secretin, and ComE (ODQ38752) (Yu et al. 2016; Klima et al. 2016).

Five GIs were predicted in the genome of M. haemolytica A2 NIVEDI/MH/1 by the Island Viewer4, which codes for 100 genes. Among these, 67 were annotated as functional genes and 33 for hypothetical proteins (Fig. 3). Comparison of M. haemolytica A2 NIVEDI/MH/1 with other M. haemolytica genome sequences using Panseq revealed that the GI3, GI4, and GI5 are unique to NIVEDI/MH/1. Among identified GIs, three regions showed higher similarity with Bibersteinia trehalosi, a bovine respiratory pathogen belonging to the Pasteurellaceae family (Blackall et al. 2007). The other two regions showed higher similarity with Actinobacillus suis, a swine pathogen (Michael J. Yaeger 1996) and Ornithobacterium rhinotracheale, a poultry respiratory pathogen (Thachil et al. 2009) (Table 3). Functional genes identified in the GIs include the plasmid stabilization system protein (ODQ39236), transposase (ODQ39232, ODQ38617), SprT family protein (ODQ38548), toxin MazF (ODQ38536), and various transcriptional regulators. It was found that GI2 has a gene for mercuric reductase (ODQ39227) needed for the mercury resistance by reduction of mercuric to elemental Hg (Freedman et al. 2012). Also, multi-copper polyphenol oxidoreductase (ODQ37742) (Guo et al. 2016), tellurite resistance methyltransferase TehB (ODQ38795), arsenate reductase (ODQ37201) (Jackson and Dugas 2003), and multidrug transporter (ODQ39218) genes were present in the GIs. Presence of these resistance genes enhances adaptability and helps the bacteria to survive in a changing hostile environment during infection. Thus, the horizontal transfer of genes among the members of the polymicrobial community might be responsible for the continuous exchange of the gene pool among the respiratory tract bacteria (Schroeder et al. 2017). HGT might lead to the development of virulence (Ho Sui et al. 2009), antibiotic resistance, and metal resistance (Hall 2010; Gilmore et al. 2013; Pagano et al. 2016).

Fig. 3

Five GIs predicted in M. haemolytica A2 NIVEDI/MH/1 using IslandViewer. The black circle represents the total base pair, the inner circle represents the GC content, and the curved line indicates the total number of contigs. Orange, blue, and red represent the predicted GIs as per SIGI-HMM, IslandPath-DIMOB, and predicted by at least one method (unstated) within the island viewer tool, respectively

Table 3 Identification of genomic islands (GIs) in M. haemolytica A2 NIVEDI/MH/1

The genome of M. haemolytica A2 NIVEDI/MH/1 also contains CRISPR-Cas array. The CRISPR-Cas is the adaptive immune system of the bacteria against invading phage and plasmid, which comprises of conserved direct repeat (DR) sequences, variable spacer sequences derived from the phage DNA, and Cas proteins (Koonin and Makarova 2009; Horvath and Barrangou 2010; Deveau et al. 2010). The CRISPR-Cas system provides immunity to bacteria against these external elements via the three-stage process of adaptation, expression, and interference (Karginov and Hannon 2010; Makarova et al. 2011). Besides the adaptive immunity against phages, the CRISPR-Cas system is also associated with gene expression and regulation of cellular processes like biofilm formation, spore formation, replicon maintenance and segregation, and DNA repair-recombination (Szczepankowska 2012). The CRISPR-Cas system is divided into three types and further into a few subtypes. Cas proteins adjacent to CRISPR sequences have the characteristic of the nucleases, helicase, polymerase, and various RNA-binding proteins (Jansen et al. 2002). Earlier, the presence of type I-C/Dvulg CRISPR-Cas system in serotypes A1, A2, and A6 of M. haemolytica were reported (Klima et al. 2016). With the use of CRISPRfinder, we identified the presence of two confirmed and two questionable CRISPR-Cas arrays in the genome of M. haemolytica NIVEDI/MH/1. Contig 4 and contig 7 of the genome contained type I-C and type II-C of the confirmed CRISPR-Cas array system, respectively. CRISPRone predicted the presence of various Cas proteins in CRISPR-Cas array. Type I-C CRISPR-Cas array of M. haemolytica NIVEDI/MH/1 consisted of Cas1, Cas2, Cas3HD, Cas4, Cas5, Cas7b, and Cas8c proteins and 23 spacers of phage DNA with 5 different direct repeats. Type II-C CRISPR-Cas array consisted of four spacers with direct repeats and Cas1, Cas2, and Cas9 proteins. Cas9 is the characteristic of the type II system and is a large protein needed for the crRNA synthesis and cleavage of the target DNA (Makarova et al. 2011). A schematic representation of the CRISPR loci along with spacers is shown in Fig. 4 (Supplementary File S4).

Fig. 4

Schematic representation of the CRISPR-Cas system in the M. haemolytica NIVEDI/MH/1 serotype A2. (a) CRISPR-Cas loci which include CRISPR array of direct repeats, a spacer sequence, and a series of Cas proteins. (b) Sequence of seven unique direct repeats (DR) and 27 spacers present in the CRISPR-Cas system present in two confirmed CRISPR array as shown in (a)

In summary, the whole genome of M. haemolytica serotype A2 isolated from a pneumonic sheep in India was sequenced. The phylogeny analysis showed that the M. haemolytica strain NIVEDI/MH/1 is closely related to other pathogenic M. haemolytica serotype A2 strains of bovine and ovine origin. The genome contained various virulence factor determinants, five GIs, and metal- and antibiotic resistance genes. Most of the GIs have shown sequence identity with other respiratory pathogens belonging primarily to the family Pasteurellaceae. Further, the strain contains type II-C of CRISPR-Cas array in its genome, in addition to the earlier reported type I-C CRISPR-Cas array. Overall, the presence of these elements substantiates the virulence potential of this Indian strain and emphasizes the role of HGT in the pathogenesis of an opportunistic respiratory pathogen. Further studies are needed on the various serotypes of M. haemolytica from different anatomical locations of the multiple hosts to understand the genotypic and pathogenic variants circulating in the Indian subcontinent.

Nucleotide sequence accession number

This Whole Genome Shotgun project has been deposited at DDBJ/EMBL/GenBank under the accession number MEHR00000000 with version MEHR01000000.