Comparative genomic analysis of clinical and environmental strains provides insight into the pathogenicity and evolution of Vibrio parahaemolyticus

Background Vibrio parahaemolyticus is a Gram-negative halophilic bacterium. Infections with the bacterium could become systemic and can be life-threatening to immunocompromised individuals. Genome sequences of a few clinical isolates of V. parahaemolyticus are currently available, but the genome dynamics across the species and virulence potential of environmental strains on a genome-scale have not been described before. Results Here we present genome sequences of four V. parahaemolyticus clinical strains from stool samples of patients and five environmental strains in Hong Kong. Phylogenomics analysis based on single nucleotide polymorphisms revealed a clear distinction between the clinical and environmental isolates. A new gene cluster belonging to the biofilm associated proteins of V. parahaemolyticus was found in clincial strains. In addition, a novel small genomic island frequently found among clinical isolates was reported. A few environmental strains were found harboring virulence genes and prophage elements, indicating their virulence potential. A unique biphenyl degradation pathway was also reported. A database for V. parahaemolyticus (http://kwanlab.bio.cuhk.edu.hk/vp) was constructed here as a platform to access and analyze genome sequences and annotations of the bacterium. Conclusions We have performed a comparative genomics analysis of clinical and environmental strains of V. parahaemolyticus. Our analyses could facilitate understanding of the phylogenetic diversity and niche adaptation of this bacterium. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-1135) contains supplementary material, which is available to authorized users.


Background
Foodborne diseases remain a serious public health problem worldwide. In Hong Kong, infections due to foodborne pathogens are also a common and important public health issue. Among all causative agents resulting in foodborne outbreaks, bacteria have caused more than 80% of the cases. In 2006, the Centre for Health Protection of Hong Kong revealed more than 800 local foodborne outbreaks due to bacterial causative agents, inflicting more than 3000 people [1]. Vibrio parahaemolyticus, a Gram-negative pathogenic halophilic bacterium, is the most prevalent in certain Asian areas and causes food poisoning, occasionally outbreaks, especially in the hot season [2].
V. parahaemolyticus is well known as the causative agent of the most prevalent food poisoning in Asia since the mackerel-borne outbreak in 1959 [3]. V. parahaemolyticus infections arise from the consumption of raw or undercooked seafood, typically causing gastroenteritis [4]. Infections can become systemic and can be lifethreatening to immunocompromised individuals [5].
V. parahaemolyticus strains isolated from diarrheal patients produce thermostable direct hemolysin (TDH), TDH-related hemolysin (TRH), or both, while isolates from the environment rarely contain genes encoding these proteins. The TDH is the major virulence gene of V. parahaemolyticus and present in most of the clinical Kanagawa phenonemon (KP)-positive strains. In addition to its hemolytic activity, enterotoxicity, cytotoxicity as well as cardiotoxicity have been confirmed in TDH [6]. TDH-related hemolysin (TRH) is produced by KPnegative strains, encoded by trh gene, with about 67% amino acid sequence homology and immunologically related with TDH [7]. TRH is hemolytic, heat labile and is suspected to play an unclarified role in the diarrhea caused by these KP-negative strains. After the bacteria produce TDH and lipopolysaccharides released from a dead cell, human host would show strong systemic responses and produce mucosal B-cells against these antigens. However, recent study has also revealed that V. parahaemolyticus could profoundly disturb epithelial barrier function in Caco-2 cells with other virulence factors, in the absence of TDH [8]. V. parahaemolyticus also encodes two type III secretion systems (T3SS), which are located in different chromosomes. The contribution of these two T3SS has previously been characterized in animal models [9,10]. These two T3SSs are involved in distinct pathogenic mechanisms during infection. The type III secretion system 1 (T3SS1) is required for cytotoxicity, whereas the type III secretion system 2 (T3SS2) is often responsible for enterotoxicity and intestinal fluid accumulation. The overall mechanism of pathogenesis in V. parahaemolyticus remains unclear.
V. parahaemolyticus infections are associated with multiple serotypes. Analysis of an outbreak in Calcutta in 1996 has identified a unique serotype clone, O3:K6, which is not previously isolated in that area [11]. After this report, O3:K6 isolates and its serovariants, defined to have identical genotype and molecular characteristics to those isolated in Calcutta, have been reported from foodborne outbreaks across the world. These serotypes were thought to be clonal derivatives of the O3:K6 serotype [12]. Recent multilocus sequence typing (MLST) data have further confirmed that multiple serotypes occur in a single genetic lineage [13].
Comparative genomics analysis of merely clinical strains of V. parahaemolyticus was described before [14]. However, the phylogenetic diversity and niche adaptation of V. parahaemolyticus at the species level is currently unknown. Strains isolated from environmental reservoirs could be as virulent as clinical strains in animal models [15], and the existing regulatory mechanisms of environmental V. parahaemolyticus could favor regulations of foreign virulence genes [16]. Nonetheless, the virulence potential of V. parahaemolyticus at the genome-wide level has not been described before. We sequenced the genomes of four O3:K6 and its serovariant clinical strains and five environmental strains of V. parahaemolyticus in Hong Kong. Genomes of a few environmental V. parahaemolyticus isolates have been announced recently [17][18][19][20], however, comparative genomics study has not been reported until now. Phylogenomic analysis of our local strains and other available genomes revealed a clear distinction between clinical and environmental strains.

Results and discussion
Summary of genome sequencing data As summarized in Table 1 obtained in Hong Kong, including four clinical and five environmental strains. The average sequencing coverage was between 44× to 66×. Genome sequences of these local V. parahaemolyticus strains were de novo assembled using Newbler 2.7 (Roche Diagnostics). Genome assembly of V. parahaemolyticus yielded 130 to 1205 contigs of N50 contig length ranging from 7 kb to 305 kb. The clinical strains we sequenced were new O3:K6 and its serovariant (O3:K59), and they were positive for the tdh gene but negative for the trh gene. Our environmental strains were all tdh and trh negative. We have also included in our analyses six publicly available V. parahaemolyticus genomes, which are all clinical strains ( Table 1) and two of them (AQ3810, AQ4037) pre-pandemic strains [14].

Phylogenetic relationships
Phylogenetic relationships among our V. parahaemolyticus isolates and the reference strains (Table 1) were examined using a genome-wide approach based on 169,998 single nucleotide polymorphisms (SNPs). Eighty nine percent of these SNPs were in coding regions, of which 25.8% were missense, 0.6% was nonsense, and 73.6% were synonymous. Genome-wide comparison readily resolved the relationship among the clinical and environmental isolates ( Figure 1). Clinical strains were found to be more related compared to the environmental strains. Moreover, our local clinical isolates were highly similar to RIMD 2210633 isolated from Thailand in 1996, consistent with previous results using various molecular methods [21] [22] and further confirm clonality of new O3:K6 and its serovariants.

Pan-genome analysis
Comparative analysis of nine newly sequenced genomes and six publicly available genomes of V. parahaemolyticus could help us determine the global gene repertoire of the species. This can be described by its 'pan-genome' that includes a core genome containing genes present in all strains and a dispensable genome composed of genes absent from one or more strains and genes that are unique to some strains. Here, the core and pan-genomes of V. parahaemolyticus were identified using OrthoMCL [23].
Our result shows that the number of pan-genome gene families increased with the number of genomes analysed, indicating that V. parahaemolyticus harbors an open pan-genome ( Figure 2). We have also analyzed the trend of new gene families, and we found that new genes will continue to be found with increasing of much more genomes.

Gene-content analysis of V. parahaemolyticus
By investigating the presence and absence of genes in V. parahaemolyticus, we found that most regions within their genomes were conserved ( Figure 3). We then further analyzed 24 genomic regions that are unique to V. parahaemolyticus RIMD 2210633 [24], which mainly include two T3SS, seven pathogenic islands, and f237 prophages. Comparative genomics revealed that these 24 regions accounted for most of the varying regions. Five pathogenic islands (VPal-1, VPal-3, VPal-4, VPal-5, and VPal-6) were found absent in not only the environmental strains, but also two pre-pandemic clinical strains. VPal-7 was found absent in AQ4037 and other environmental strains.
Previously, microarray-based approach failed to detect genes responsible for pathogenesis of V. parahaemolyticus [25]. Traditional features also failed to differentiate clinical strains from environmental strains [26,27]. Here, analysis of these varying regions revealed the presence of only 11 regions in the clinical group, which were absent in the environmental group ( Table 2). As expected, genes in T3SS2 were not found in all of the clinical strains, which further confirm that it is not reliable to simply use tdh, trh or genes in T3SS as virulence features [25]. However, a region located in chromosome II (984909-999901) contained genes from VPA0950 to VPA0957, which encode putative biofilm associated proteins and outer membrane proteins. A few genes including two type IV pili were previously described in biofilm formation of V. parahaemolyticus, however details of the process are still unexplored in this bacterium [28,29]. In other Vibrio species such as V. fischeri, biofilm formation was found to play an important role in  host colonization [30]. It is suggested that gene clusters found in clinical strains of V. parahaemolyticus could also provide this bacterium with advantages in host adaptation. These gene clusters have not been found previously, and we suggest that they could act as new candidate virulence markers.
A novel genomic island commonly found in clinical strains O3:K59 is a new serovariant of serotype O3:K6. It is recently reported that it could replace local serovariants and has caused pandemics in Chile and Peru [31,32].
Here we analyzed the genome of our O3:K59 strain VIP4-0407 and found a novel genomic island that is absent in V. parahaemolyticus RIMD 2210633 [33]. This genomic island was found in most of our clinical strains including AQ4037, K5030, AN-5034, Peru-466, VIP4-0439, VIP4-0395, VIP4-0407 and AQ3810, but was absent in our environmental strains. We defined this genomic island as Vibrio parahaemolyticus island-8 (VPal-8). It contained six CDSs (Table 3), one of which showing high similarity with an AraC family transcriptional regulator that caused virulence in a mouse infection model in Mycobacterium tuberculosis [34]. ArcC family transcriptional regulators could regulate T3SS genes to modulate bacterial virulence [35], We analyzed sequences in this genomic island to determine whether there is any T3SS secreted protein. Using a combination of available T3SS prediction tools including T3_MM [36], Effective T4 [37], and T3SS effector prediction [38], we found that the third protein in this genomic island showed positive results in all of the prediction tools.
We further investigated whether VPal-8 is widespread across the population of clinical strains. Using this genomic island as a probe to perform southern hybridization on an additional population of 40 clinical strains, we found that this novel genomic island can be found in 29 of the clinical strains (Additional file 1: Table S1). Therefore, we suggest that the genomic island found in this study may be related to the virulence of V. parahaemolyticus.

sRNAs and CRISPR element analysis
CRISPRs (Clustered regularly interspaced short palindromic repeats) could play important roles in the interaction of bacteria and mobile genetic elements [39]. We annotated CRISPR elements in V. parahaemolyticus using CRISPR finder [40] and found a total of six CRISPR elements (Table 4). Our environmental strains were found to harbor fewer CRISPR types, whereas at least two CRISPRs were found in each of the clinical strains.
sRNAs are a class of non-coding RNAs in bacteria and are important post-transcriptional regulators in multiple crucial biological processes such as biofilm formation, quorum sensing and virulence [41,42] . We determined the sRNA sequences in our V. parahaemolyticus strains and found that most sRNAs were highly conserved (Additional file 2: Table S2). However, a few sRNAs such as GcvB and STnc1460 were only present in a few strains, indicating that these sRNAs may mediate some strain-specific regulations.

Biphenyl degradation pathway identified in environmental strains
We annotated the subsystems of V. parahaemolyticus genomes using the RAST Server [43] (Figure 4). Five subsystems were found enriched in the environmental strains, including Metabolism of Aromatic Compounds (p = 1.12E06), Protein metabolism (p = 0.01), Iron acquisition and metabolism (p = 0.027), Phages and Prophages (p = 0.041) and Sulfur Metabolism (p = 0.046). The most enriched subsystem "Metabolism of Aromatic Compounds" contained an average of six genes in clinical strains, but an average of 10 genes in the environmental strains.
We further examined this category to investigate which pathway has contributed most of the differences. Genes involved in biphenyl degradation pathway (BphE1, BphC, Bphj2) were found in all of our environmental strains but not in the clinical strains, except the prepandemic strain AQ3810 (Table 5). The presence of the unique biphenyl degradation pathway in our environmental strains is possibly caused by the widespread of polychlorinated biphenyls in the environment [44].
Virulence potential of environmental strains V. parahaemolyticus harbors many virulence factors including TDH, TRH and two T3SSs, which are generally pathogenic. Here we report the distribution of virulence  factors and prophages in the genomes of our V. parahaemolyticus strains to reveal their pathogenic potential. We predicted prophage regions in V. parahaemolyticus using the PHAST webserver [45] and found that clinical strains were restricted to a few prophage elements. The pandemic group of clinical strains possessed a filamentous prophage f237, which could increase virulence by adhering to intestinal cells [46]. A similar prophage was also found in the AQ3047 strain. In contrast, many diverse prophage elements were found in our environmental strains (Table 6). Six new prophages were found in four of our environmental strains, which harbor genes that could enhance the fitness for survival in specific environments. We further examined the relationship between the number of prophages and CRISPR in our strains. We analyzed the spacer sequences of CRISPRs and eight unique spacers sequence were found. However, we could not find any homologous prophage regions that could match with these spacer sequences in the respective strains. This could be explained by the possibility that the CRISPRs of Vibrio   parahaemolyticus could inhibit the prophage insertion to the host genome. Indeed, an inverse correlation of CRISPRs with the number of prophage was also described in Streptococcus pyogenes [47]. toxR gene, located next to toxS gene, involves in the regulation of many virulence-associated genes of V. parahaemolyticus [48]. Variation in the toxR sequence can be used to differentiate phylogenetically distinct clusters in V. parahaemolyticus. Seven base positions were previously identified to distinguish the O3:K6 isolates before 1995 from the new O3:K6 clones within a 1346-bp region [49]. We further investigated the variations between these base positions within our V. parahaemolyticus isolates. In all of the seven new O3:K6 clones, the bases were perfectly mapped to the reference ones (Table 7). However, another environmental strain (VIP4-0443) also had the same sequence bases with the new O3:K6 clones, suggesting that it may also harbor virulence potential. We then further examined whether there are any other virulence proteins in this strain. A total of 264 candidate virulence genes were found (Additional file 3: Table S3), which further confirm that this environmental strain could have virulence potential.
Construction of a genome sequence database for V. parahaemolyticus A V. parahaemolyticus genome database was constructed here based on the Ensembl genome annotation system using Perl scripts ( Figure 5). Genome DNA sequences can be input into the platform via a user-friendly web-based interface. Our web-based database allows comparative analysis and data mining using available V. parahaemolyticus genome data. A suite of useful computational tools is now available for data analysis such as MLST and detection of genetic variations. The web-based database is now launched as a publicly accessible domain (http://kwanlab. bio.cuhk.edu.hk/vp).

Conclusions
This is the first study on the genome dynamics of V. parahaemolyticus at the species level. We have three main conclusions. First, by comparing the genome sequences, we found that our clinical strains were phylogenetically distinct from the environmental strains. We discovered a new gene cluster belonging to the biofilm associated proteins of V. parahaemolyticus. We suggest that this novel gene cluster is a new virulence marker differentiating clinical from environmental strains. We also found a novel genomic island (VPal-8) that frequently distributed in the clincial strains. Second, by analyzing the virulence features and prophage elements of our environmental strains on a genome-wide scale, we discovered a new type of toxRS in one of the environmental strains and a diverse array of prophage elements in the environmental strains. These results indicate that the environmental strains studied could have virulence potential. Third, we set up an Ensemblbased genome database for V. parahaemolyticus to provide a unified and user-friendly platform to access all currently available genome sequences and annotations of this bacterium.

Ethics statement and isolate collection
Nine isolates of V. parahaemolyticus were selected from the bacterial archive maintained at the microbiology laboratory of the Department of Health of Hong Kong and were characterized for genome sequencing. The study was approved by the Joint CUHK-NTEC Clinical Research Ethical Committee at the Prince of Wales Hospital in Hong Kong. Written informed consent was obtained from all the studied subjects for sample collection and subsequent analysis.

DNA sequencing
Whole-genome shotgun sequencing was performed using a 454 Genome Sequencer (GS) FLX-Titanium (Roche Diagnostics, US). DNA of each selected isolate was extracted from culture and subjected to sequencing library preparation according to the manufacturer's recommended protocols. Libraries of processed genomic DNA fragments immobilized on DNA capture beads were individually sequenced on a PicoTiterPlate device. Bases sequenced and the corresponding quality values were called and delivered in a standard format by GS-FLX-Titanium system for downstream bioinformatic analyses. The remaining genome gaps were filled by the GS-FLXprovided paired-end approach and/or primer walking. The primer walking approach involved design of primers flanking the gaps and PCR to amplify the DNA fragments covering the gaps. PCR products were then directly sequenced using the ABI DNA Sequencer (Applied Biosystems, US).

Sequence assembly and annotation
Raw sequence reads were first filtered to remove lowquality reads. Sequence assembly was performed on the     A  A  T  T  T  A  T   Peru 466  Clinical  A  A  T  T  T  A  T   K5030  Clinical  A  A  T  T  T  A  T   AN-5034  Clinical  A  A  T  T  T  A  T   VIP4-0439  Clinical  A  A  T  T  T  A  T   VIP4-0395 Clinical remaining clean reads de novo using the GS Newbler 2.7. The genome sequences were annotated using the RAST webserver. Contigs were reordered and genome comparisons were carried out using mauve program [50]. BRIG was used for genome alignment visualization [51]. Prophage elements of V. parahaemolyticus were predicted using PHAST [45].

Phylogenetic tree construction
SNPs among our isolates and several foreign strains were identified using GS Reference Mapper (Roche Diagnostics). SNPs detected among all strains were concatenated together and a neighbor-joining phylogenetic tree was constructed using MEGA 6 [52] with 500 bootstrap pseudoreplicates. Vibrio harveyi ATCC BAA-1116 was used as outgroup for defining the root.

Small RNAs and CRISPRs
A list of small RNA sequences was retrieved from the BSRD database [42]. sRNA homologs in our V. parahaemolyticus isolates were then identified using BLASTn. The e-value of BLAST was set to 1e-5. Annotation of CRISPR elements was done using CRISPRfinder [40].

Experimental validation of genomic island
Two sets of primers were designed for PCR-based detection of genomic island in a population of clinical V. parahaemolyticus. These clinical strains were previously used for microarray-based analysis [53].