Complete genome sequence and comparative genomics of the golden pompano (Trachinotus ovatus) pathogen, Vibrio harveyi strain QT520

Ethics statement

All protocols for experiments involving live animals conducted in this study were approved by the Animal Experimental Inspection of Laboratory Animal Centre, Hainan University (Haikou, China), and the approval number is 20160709. We conscientiously abided by the ethical principles of animal welfare and the rule of laboratory animal center.

Bacterial strain isolation

V. harveyi QT520 was initially isolated in Qiaotou, China from T. ovatus showing clinical signs of body rot. To obtain a pure culture, a single colony of strain QT520 was selected and transferred to fresh TCBS medium at least three times. Then, the clone was cultivated at 30 °C in 2216E liquid medium for 16 h. The cell morphology of QT520 and contamination of the culture with other microorganisms were assessed by light microscopy (BA410; Motic, Xiamen, China). The strain was archived at Hainan University and Hainan Academy of Ocean and Fisheries Sciences under accession number QT520.

Genomic DNA extraction and whole genome sequencing

Genomic DNA was extracted using a Rapid Bacterial Genomic DNA Isolation Kit (Sangon Biotech, Shanghai, China) according to the manufacturer’s instructions. The concentration and purity of the extracted DNA were determined using a DNA Qubit 2.0 (Invitrogen, Carlsbad, CA, USA). The extracted genomic DNA was used to construct next-generation libraries with an insert size of 500 bp using a NEBNext^® Ultra™ DNA Library Prep Kit for Illumina^® and PacBio RS II libraries with an insert size of 10 kb.

Genome sequencing was conducted by Sangon Biotech Co., Ltd. (Shanghai, China) using the next-generation Illumina MiSeq (Illumina, San Diego, CA, USA) and the third-generation PacBio RS II sequencing platform (Pacific Biosciences, Menlo Park, CA, USA), respectively. The raw sequences obtained from the third-generation PacBio RS II sequencing platform were assembled using Canu (Koren et al., 2016). Gapcloser and GapFiller were used to close the gaps with next-generation sequence reads where possible after assembly (Boetzer & Pirovano, 2012) and the sequence variants were detected and assembled by PrInSeS-G to obtain high-quality, whole-genome sequences (Massouras et al., 2010).

Genome analysis

Prokka was used to predict and annotate prokaryotic ORFs, tRNAs and rRNAs (Seemann, 2014). Clusters of Orthologous Groups (COG) analysis was performed using RPS BLAST. Antibiotic resistance genes and genes encoding virulence factors were identified through BLAST searches of the Comprehensive Antibiotic Resistance Database (CARD) and the Virulence Factors of Pathogenic Bacteria Database (VFDB), respectively. Circular genome maps were generated using the CGView Server (Grant & Stothard, 2008) based on the information generated by the genome annotation.

Phylogenetic analyses

Four phylogenetic trees were constructed based on average nucleotide identity (ANI), core genes, multilocus sequence typing (MLST), and the 16S rRNA gene, respectively. The ANI values were calculated using Jspecies (Goris et al., 2007) and phylogenetic tree analysis of 12 available complete genome sequences of Vibrio sp. strains (Table S1) based on ANI values was performed using the R package Ape with neighbor-joining methods (Paradis, Claude & Strimmer, 2004). The core genes of the 12 genomes were obtained using cd-hit 4.6.1, the corresponding protein sequences were aligned using MUSCLE 3.8.31 and the phylogenetic trees were generated using Treebest 1.9.2 (Caputo et al., 2015; Mbengue et al., 2016). Four different housekeeping genes were used to construct the MLST phylogenetic tree: toxR, vhhA, ompK and hsp60. The 16S rRNA gene sequences and the corresponding amino acid sequences inferred from the four genes for MLST were respectively aligned with those of the most closely related species using the multiple alignment program ClustalW. Phylogenetic relationships between strain QT520 and closely related species were determined using the MEGA6 software. Phylogenetic trees were generated using the neighbor-joining (NJ) method with 1,000 randomly selected bootstrap replicates.

Comparative genomics

Genomes of ATCC 33843 (392 (MAV)), ATCC 43516, and QT520 were compared through Mauve using the default parameters (Darling et al., 2004; Mbengue et al., 2016). The locations of different genes in the gene clusters were visualized using SVG 1.1, and ImageMagick 6.5.4–7 was used to convert SVG format to PNG.

Determination of Median Lethal Dose (LD50)

The isolated QT520 strain was inoculated into 2216E liquid medium and incubated on a shaker (180 rpm) at 30 °C overnight, after which a culture containing 1.47 × 10⁸ CFU/mL obtained. During challenge, the bacterial suspension was serially diluted 10-fold using Stroke-Physiological Saline Solution and was used to infect T. ovatus (10 fish/group) by intraperitoneal injection with 1.47 × 10⁸ CFU/mL, 1.47 × 10⁷ CFU/mL, 1.47 × 10⁶ CFU/mL or 1.47 × 10⁵ CFU/mL. The control group was injected with Stroke-Physiological Saline Solution. Fish deaths were recorded for 2 weeks and the Reed-Muench method was used to calculate the median lethal dosage.

General features of V. harveyi QT520

Through the use of a combination of sequencing approaches using the MiSeq and PacBio platforms, we were able to assemble five contigs with 227.76 × coverage. The genome sequence of strain QT520 was deposited in the DDBJ/EMBL/GenBank under the accession numbers CP018680–CP018684. The project information, according to the minimum information about a genome sequence (MIGS) recommendation (Field et al., 2008), is shown in Table S2. The strain QT520, which was isolated from cultured T. ovatus exhibiting clinical signs of body rot, was Gram-negative and aerobic. After negative staining it was observed to be an oval, rod-shaped cell approximately 1.6–2.2 µm in length and 1.0–1.1 µm in width with a single polar flagellum (Fig. 1).

Genomic information

The QT520 genome was a total of 6070846 bp, consisting of two circular chromosomes and three circular plasmids (chromosome I, 3,560,044 bp; chromosome II, 2,260,627 bp; plasmid p1, 113,574 bp; plasmid p2, 76,744 bp; plasmid p3, 59,857 bp) (Fig. 2). The genome had an overall G  + C content of 45% (Table 1). A total of 3,377 open reading frames (ORFs) were detected in Chromosome I, including 3,195 coding sequences (CDSs), 118 tRNA genes, 34 rRNA genes, 29 miscRNA and 1 tmRNA; 2,425 CDSs (75.90% of 3,195 CDSs) were annotated as functional genes, and 770 CDSs (24.10%) were annotated as hypothetical or uncharacterized genes. Chromosome II contained 2,061 ORFs, 2,030 CDSs, 16 tRNA genes, three rRNA genes and 12 miscRNA; 1,476 CDSs (72.71%) were annotated as functional genes, and 554 CDSs (27.29%) were annotated as hypothetical or uncharacterized genes. The chromosome maps (Fig. 3) display 3,657 ORFs (64.15% of the total number of predicted ORFs) that encode known functional proteins, and 423 ORFs having no known function in our COG functional categorization. Among the functionally predicted ORFs, 2,593 ORFs (63.55% of the COG-assigned ORFs) belonged to nine major COG functional categories: 527 ORFs in category R (General function prediction only), 329 ORFs in category K (Transcription), 287 ORFs in category E (Amino acid transport and metabolism), 276 ORFs in category T (Signal transduction mechanisms), 254 ORFs in category G (Carbohydrate transport and metabolism), 241 ORFs in category C (Energy production and conversion), 243 ORFs in category P (Inorganic ion transport and metabolism), 236 ORFs in category M (Cell wall/membrane/envelope biogenesis) and 200 ORFs in category J (Translation, ribosomal structure and biogenesis).

Virulence factors of V. harveyi QT520

Virulence factors of strain QT520 are listed in Table S3. The results of the QT520 genome analysis revealed a number of common virulence determinants, such as ACF, IlpA, OmpU, Flagella, T6SS, TTSS, Cya, Hemolysin and MARTX. Factors responsible for iron acquisition, including IutA, Chu, Enterobactin, FeoAB and Mycobactin were also present. In addition, some adhesion proteins (Capsule/CapsuleI, Los, IlpA, MAM7 and OmpU), antiphagocytosis proteins (Alginate), efflux pumps, endotoxins, and stress proteins were also found in strain QT520.

Phylogenetic analysis

Four comparative phylogenetic tree analysis methods were used to identify closely related strains of QT520: ANI, core genes, MLST, and 16S rRNA gene sequencing. For phylogenetic tree analysis using ANI, 11 complete genome sequences available in the NCBI GenBank database, as well as that of QT520, were obtained and the ANI values were calculated (Table S4). A phylogenetic tree based on ANI values was constructed (Fig. 4A). The tree identified a subgroup that contained 1114GL, LMB29, ATCC BAA-1116, ATCC 33843 (392 (MAV)), ATCC 43516 and QT520, which was designated Group I. Within this group, the genome of QT520 was most closely related to those of ATCC 33843 (392 (MAV)) and ATCC 43516 (98.49 and 98.59% of ANI values, respectively).

An additional phylogenetic tree of the 12 complete genome sequences described above was constructed using the core genes method (Fig. 4B). The results showed that ATCC 33843 (392 (MAV)) and ATCC 43516 were the closest in Group I and that QT520 was closely related to both strains, supporting that the ANI and core genes methods could distinguish them even at the strain level.

Another phylogenetic tree was generated using MLST with four different housekeeping genes: toxR, vhhA, ompK and hsp60. The DNA sequences of the four MLST genes in the Group I strains were compared and aligned to generate a phylogenetic tree (Fig. 4C). The results showed that the genomes of QT520, ATCC 33843 (392 (MAV)) and ATCC 43516 were again the closest relatives and those of 1114GL, LMB29 and ATCCBAA-1116 were the closest in Group I. However, comparative phylogenetic tree analysis using 16S rRNA sequences of the Group I strains revealed that strains QT520, 1114GL, ATCC 33843 392 MAV and ATCC 43516 were the closest relatives, being classified into one phylogenetic branch (Fig. 4D).

Comparative genomics of V. harveyi ATCC 3843 (392 (MAV)) and ATCC 43516

Comparative genomic analyses were conducted to compare V. harveyi QT520 with its most closely related strains, V. harveyi ATCC 3843 (392 (MAV)) and ATCC 43516 (Fig. 5A). The results showed that the number of core genes among the three strains is 4,463. The number of specific genes in QT520, ATCC 3843, and ATCC 43516 are 458, 277 and 448, respectively.

Two genes (APP0488.1 and APP06774.1) identified as core genes encode a multi-antimicrobial extrusion efflux family protein and a xanthine-guanine phosphoribosyltransferase, respectively. The results suggested that the three V. harvey i strains may have multiple antibiotic resistances to tigecycline, tetracycline, streptomycin, kanamycin, ciprofloxacin and norfloxacin (Kim et al., 2003; Kuroda & Tsuchiya, 2009; Teo, Tan & Poh, 2002).

Virulence factor genes, including APP05600.1, APP06379.1, APP06380.1, APP06382.1, APP06383.1, APP06387.1, APP06388.1, APP07902.1, APP09130.1, APP09132.1, APP09134.1, APP09207.1 and APP09259.1 were only observed in strain QT520 based on the comparative analysis results (Table S3 and Fig. 5A). These specific genes encode Type IV pili, O-antigen, LPS, Capsule, Cya, hemolysin, MARTX and TTSS (SPI-1), respectively. These genes are well known to play important roles in attachment, expression of other virulence factors, preventing phagocytosis, anti-inflammatory effects, hemolytic activities and a cytotoxicity (Ruwandeepika et al., 2012; Zhang & Austin, 2000), suggesting that they may be responsible for the virulence associated with the QT520 strain.

Genes associated with type I secretion systems (T1SS) and type VI secretion systems (T6SS) were detected in QT520, but not in ATCC 33843 (392 (MAV)) and ATCC 43516. The numbers of type II, III and IV secretion systems (T2SS, T3SS, and T4SS) related genes in QT520 was 28, 32, and 21, respectively; the numbers of T2SS, T3SS, and T4SS related genes in ATCC 33843 (392 (MAV)) was 7, 12, and 4, respectively. The numbers of genes T2SS, T3SS, and T4SS related genes in ATCC 43516 was 7, 14, and 4, respectively. Greater numbers of genes associated with the T1SS, T2SS, T3SS, T4SS and T6SS were detected in the genome sequence of QT520 than in strains ATCC 33843 (392 (MAV)) and ATCC 43516 (Table 2 and Fig. 5). In addition, T1SS were detected in regions of chromosome I (APP04465.1), chromosome II (APP08059.1 and APP08730.1) and plasmid p1 (APP09130.1 and APP09132.1) of QT520. T2SS and T3SS were detected in the regions APP04623.1–APP08089.1 and APP04087.1–APP08428.1 in chromosomes I and II, respectively. Nine T4SS-encoding genes were detected in the insertion region in QT520 compared to ATCC 3843(392(MAV)) and ATCC 43516 (Fig. 5B and Fig. 5C), and other T4SS-encoding genes were primarily distributed in regions of chromosome I (APP05190.1–APP05200.1) and chromosome II (APP07457.1 to APP08116.1). T6SS-related genes were detected in regions of chromosome I (APP04433.1 to APP04446.1) and chromosome II (APP07649.1 to APP07972.1).

Three plasmids were detected in V. harveyi strain QT520 that were not observed in strains ATCC 3843 (392 (MAV)) and ATCC 43516. In particular, three important virulence factor genes (cyaB, hlyB and rtxA) were present in plasmid p1.

LD₅₀

Artificial infection results showed that the LD₅₀ of QT520 was 2.5 × 10⁵ bacteria per fish, suggesting that this bacterium is strongly pathogenicity to T. ovatus (Table 3).