ABSTRACT
Vibrio parahaemolyticus is a bacterial pathogen that can cause illness after the consumption or handling of contaminated seafood. The primary virulence factors associated with V. parahaemolyticus illness are thermostable direct hemolysin (TDH) and Tdh-related hemolysin (TRH). However, clinical strains lacking tdh and trh have recently been isolated, and these clinical isolates are poorly understood. To help understand the emergence of clinical tdh- and trh-negative isolates, a genomic approach was used to comprehensively compare 4 clinical tdh- and trh-negative isolates with 16 environmental tdh- and trh-negative isolates and 34 clinical isolates positive for tdh or trh, or both, with the objective of identifying genomic features that are unique to clinical tdh- and trh-negative isolates. The prevalence of pathogenicity islands (PAIs) common to clinical isolates was thoroughly examined in each of the clinical tdh- and trh-negative isolates. The tdh PAI was not present in any clinical or environmental tdh- and trh-negative isolates. The trh PAI was not present in any environmental isolates; however, in clinical tdh- and trh-negative isolate 10-4238, the majority of the trh PAI including a partial trh1 gene was present, which resulted in reclassification of this isolate as a tdh-negative and trh-positive isolate. In the other clinical tdh- and trh-negative isolates, neither the trh gene nor the trh PAI was present. We identified 862 genes in clinical tdh- and trh-negative isolates but not in environmental tdh- and trh-negative isolates. Many of these genes are highly homologous to genes found in common enteric bacteria and included genes encoding a number of chemotaxis proteins and a novel putative type VI secretion system (T6SS) effector and immunity protein (T6SS1). The availability of genome sequences from clinical V. parahaemolyticus tdh- and trh-negative isolates and the comparative analysis may help provide an understanding of how this pathotype is able to survive in vivo during clinical illness.
INTRODUCTION
Vibrio parahaemolyticus is a Gram-negative, halophilic bacterium that is ubiquitous in marine and estuarine environments and is often found colonizing shellfish or shrimp. While most strains are nonpathogenic, many have acquired virulence factors that result in illness when individuals are exposed to V. parahaemolyticus strains carrying these virulence factors (1). V. parahaemolyticus is recognized to be a leading cause of foodborne illness worldwide and is transmitted via the handling and consumption of raw or undercooked contaminated seafood (1). Infections occur both sporadically and in very large outbreaks. The most common manifestation of V. parahaemolyticus infection is acute, watery diarrhea accompanied by abdominal pain and nausea, although symptoms can also be severe and include a dysentery-like illness or septicemia (2). Since most cases of illness caused by V. parahaemolyticus are self-limiting, rates of infection are probably underestimated due to underreporting.
Clinical V. parahaemolyticus isolates generally have at least one of two major toxigenic virulence factors, thermostable direct hemolysin (TDH) (3) and TDH-related hemolysin (TRH) (4). TDH has hemolytic activity on a blood-containing medium, Wagatsuma agar, and the process is referred to as the Kanagawa phenomenon (KP) (3). During infection, TDH is involved in cytotoxicity and hemolytic activity, and on the basis of the sequence similarity between TDH and TRH, TRH is believed to act similarly (5–7). The presence of tdh and/or trh is common in pathogenic isolates but relatively rare in environmental strains; therefore, the presence of these genes is used to assess the virulence potential of V. parahaemolyticus isolates (8, 9). In addition to tdh and trh, whole-genome sequencing (WGS) of V. parahaemolyticus led to the identification of two nonredundant type III secretion system (T3SS) gene clusters, dubbed T3SS1 and T3SS2, on chromosome 1 and chromosome 2, respectively, which are also involved in virulence (5, 10). T3SSs are needle-like apparatuses that inject bacterial effector proteins, such as toxins or hemolysins, directly through the membrane and into the cytoplasm of eukaryotic cells (11). On the basis of its G+C content and its high degree of sequence identity with the T3SSs of other Vibrio species, T3SS1 appears to have been ancestrally acquired and is present in all V. parahaemolyticus isolates, even nonpathogenic strains (10). Once it is in vivo, T3SS1 appears to inject effectors, such as VopQ, VopR, VopS, and VPA0450, directly into eukaryotic cells, resulting in cytotoxicity (5, 10). T3SS2 derives from two separate lineages: T3SS2α is typically found on a pathogenicity island (PAI) with the tdh gene, and T3SS2β is found with the trh gene. T3SS2 appears to inject VopA, VopC, VopL, and VopT into eukaryotic cells, resulting in both cytotoxicity and enterotoxicity (5, 12–14). T3SS2 has a lower G+C content than the genomic average, which is indicative of a recent acquisition of this region through lateral gene transfer (15).
Some V. parahaemolyticus strains also have two type VI secretion systems (T6SSs) (16). The T6SS was recently defined functionally in V. cholerae and is structurally similar to a contractile phage tail, but in the reverse orientation, that is fully assembled inside the bacterial cell and injects effectors directly into the recipient cell (17, 18). It is composed of 13 essential genes and a variable number of nonessential genes, including various effectors (19). In V. parahaemolyticus, T6SS2 is found in all strains, while T6SS1 is mostly associated with clinical isolates (20) and may play a role in virulence (21), although this has not yet been demonstrated conclusively. The T6SS1 also appears to have antibacterial activity in the environment, which may give isolates containing this system a competitive advantage (16).
Recent studies have reported the isolation of strains that lack both tdh and trh (tdh- and trh-negative strains) from clinical samples (22–27). However, the ability of tdh- and trh-negative strains to independently cause clinical illness is still controversial. For example, coinfection with multiple V. parahaemolyticus strains is known to occur, and if at least one infecting strain carries tdh or trh, it is possible that nonpathogenic tdh- and trh-negative strains could be isolated from a sick individual without having a direct involvement in illness (24). The coinfection model is supported by the finding that a single seafood sample is often contaminated by several different V. parahaemolyticus strains, some of which appear to be nonpathogenic (24). Alternatively, there are several lines of evidence that support the opinion that some tdh- and trh-negative isolates are able to induce clinical infection. During a coinfection study, three sick patients produced 30 tdh- and trh-negative isolates, and despite multiple culturing attempts, no other enteric pathogen or tdh- and trh-positive V. parahaemolyticus strain could be isolated from these patients (24). However, regardless of their independent pathogenicity, clinical tdh- and trh-negative isolates have a demonstrated ability to survive in vivo during illness, and we do not understand if the role of these isolates in human illness is as a causative agent, an innocent bystander, or an active participant in a multistrain infection. Therefore, this investigation was undertaken to better understand clinical tdh- and trh-negative isolates and how they compare on a genomic level to traditional pathogenic V. parahaemolyticus isolates and environmental tdh- and trh-negative isolates. Here we present a thorough comparative genomic analysis of multiple clinical tdh- and trh-negative isolates. This comprehensive approach has provided several insights into the pangenomics of clinical V. parahaemolyticus isolates and led to the identification of a novel putative T6SS effector.
MATERIALS AND METHODS
Genome sequencing, assembly, and annotation.Each of the clinical isolates sequenced in this study were Canadian clinical strains originating from provincial public health laboratories and submitted to the National Microbiology Laboratory (Public Health Agency of Canada), the British Columbia Centre for Disease Control (BCCDC), or the Bureau of Microbial Hazards (BMH) (Health Canada). All isolates were routinely propagated on TSA-2N agar (Difco BD, NJ, USA). The environmental isolates were each isolated from seafood by the Vibrio Reference Laboratory in Canada. DNA for whole-genome sequencing was extracted using a Maxwell 16 SEV cell DNA purification kit (Promega, Madison, WI). The short-read sequence data were generated by preparing a paired-end library with a Nextera XT DNA sample preparation kit (Illumina, San Diego, CA) and sequencing the library on a MiSeq benchtop sequencer (Illumina) for 500 or 600 cycles. Previously, Banerjee et al. (26) performed PCR amplification of the tdh and trh genes, and this allowed us to select isolates for WGS and group isolates according to clinical pathotype (tdh and trh positive, tdh negative and trh positive, tdh positive and trh negative, or tdh and trh negative). This grouping scheme was used during the analysis in our study.
Genomes were assembled de novo using the SPAdes (version 3.1.1) program, and quality was assessed using QUAST, as outlined by Ronholm et al. (Table 1) (47). De novo assembly resulted in the PAIs of interest, including the tdh PAI and the trh PAI, being split between contigs. For analysis of pathogenicity islands, raw MiSeq reads were reassembled by reference-guided assembly using the Burrows-Wheeler aligner (BWA; version 07.05a) (29) and then visualized and inspected using the Tablet (version 1.14.10.20) program (30).
Sequencing and annotation results for V. parahaemolyticus strains analyzed during this studya
Annotation of the function of the protein-coding sequences was performed using the Basic Local Alignment Search Tool (BLAST), which was used to compare the sequences obtained against the sequences in the COG (Clusters of Orthologous Groups) of proteins databases and the NCBI nr protein database.
Comparative genomics.An in-house Perl script was used to identify all orthologous, accessory, and unique genes between Vibrio genomes (the source code can be found at https://github.com/bfssi-nicholas-petronella/HARDCore.git.) Using the BLAST program, genes having 60% sequence identity over 80% of their length were considered orthologs (31, 32). The set of orthologous genes shared by all genomes was defined as the core genome, total genes identified within all genomes were defined as the pangenome, accessory genomes were defined as the set of genes possessed by a subset of genomes (pathotypes), and unique genomes were the subset of genes possessed only by a single strain's genome. To construct whole-genome rarefaction curves, an in-house script was also used in an iterative manner to obtain the number of genes contained in the pangenome relative to the number of genomes analyzed. For example, for the data points for the 16 genomes, the size of the pangenome was calculated for 16 randomly selected genomes (from the 38 possibilities). This was repeated 1,000 times using 16 genomes randomly selected from 38 genomes, and the average was taken. The purpose was to determine if V. parahaemolyticus had an open or closed pangenome and verify that the genomes from the number of strains selected for the experiment could adequately represent the core genome.
Phylogenetics.To construct phylogenetic trees for the clinical isolates, the largest representative sequence of each core gene was retrieved from the pangenome. The homologue for each representative sequence was then retrieved from the whole genome of each strain. The Prokka software tool was used to generate the corresponding amino acid sequence (33). Trees were constructed from the core genome amino acid sequence. Each gene was aligned using the Muscle (version 3.8.31) program with default parameters (34), and then RAxML software was used to construct a maximum likelihood phylogeny from concatenated alignments with 100 bootstrap values (35).
To construct phylogenetic trees for each isolate used in the study, the multilocus sequence typing (MLST) sequences were retrieved from the WGS. Trees were constructed from the MLST alleles. Each gene was aligned using the Muscle (version 3.8.31) program with default parameters (34), and then RAxML was used to construct a maximum likelihood phylogeny from concatenated alignments with 1,000 bootstrap values (35).
Genetic islands and secretion systems.The presence and absence of the tdh and trh pathogenicity islands and various secretion systems were investigated using BLASTn. The EasyFig (version 2.1) program was used for visualization of genomic islands (36). For larger genomic islands (tdh PAI and trh PAI), a Burrows-Wheeler transform reference-guided assembly was used on raw fastq data with the canonical sequence of the genomic island (29). The assembly was inspected for accuracy and validity using Tablet, prior to visualization (30). To possibly identify novel pathogenicity islands in the tdh- and trh-negative pathotype, the completed results of an IslandViewer analysis on strain CDC_K4557 were downloaded from the IslandViewer3 website (37). BLASTn was used to search each genome included in this study (Table 1) for each island that was identified.
Nucleotide sequence accession numbers.All nucleotide sequence data referred to in this article have been deposited in the DDBJ/EMBL/GenBank database under the various BioProject accession numbers provided in Table 1. Additional important data are included in the supplemental material.
RESULTS AND DISCUSSION
Type III secretion systems and tdh and trh pathogenicity islands.Traditionally, a pathogenic V. parahaemolyticus strain has been defined by the presence of tdh or trh, or both (5). These two virulence factors generally occur near a T3SS. Two T3SSs have been reported in V. parahaemolyticus, and two variants of T3SS2 (T3SS2α and T3SS2β) have been described (12). T3SS2α is associated with tdh, while T3SS2β is typically found with trh (14), though exceptions exist (27). T3SS1 is composed of 42 genes (VPA1656 to VAP1696 and VPA0450 in V. parahaemolyticus RIMD2210633 [15]). T3SS1 was present in each of the clinical and environmental isolates in this investigation. The tdh PAI (also known as VPAI-7) is composed of 87 coding sequences (in V. parahaemolyticus RIMD2210633) and includes tdhA, tdhS, and the T3SS2α genes (38). In our current work, due to multiple homologous areas in the tdh PAI and trh PAI areas, the use of de novo assembly led to this large PAI being split between multiple contigs. Therefore, to analyze this region properly, a reference-guided assembly based on the tdh PAI sequence (VPA1310 to VPA1396) from V. parahaemolyticus RIMD2210663 was used (39). Using this method, a complete tdh PAI including T3SS2α was consistently observed in all clinical tdh-positive and trh-negative isolates (Fig. 1A). Clinical tdh- and trh-positive isolates contained homologues to some coding sequences typically found in this island (VPA1310, VPA1311, VPA1313 to VPA1318, VPA1320, VPA1321, VPA1329, VPA1342, VPA1347, and VPA1382 to VPA1393), and their presence was consistent among all clinical tdh- and trh-positive isolates (see Fig. 3A). Isolates categorized as tdh negative and trh positive or tdh and trh negative did not contain a PAI with homology to tdh PAI genes (Fig. 1A).
Pathogenicity islands and T3SS2. tdh PAI (A) and trh PAI (B) are the two canonical genomic islands present in clinical isolates of V. parahaemolyticus. The presence or absence of these islands is strongly correlated with the presence of the tdh and trh genes. In three of the four tdh- and trh-negative isolates, these genes are entirely absent. (C) In the fourth tdh- and trh-negative isolate, 10-4328, a trh PAI including a partial trh1 gene is mostly present. Abbreviations: neg, negative; pos, positive.
The trh PAI, composed of 81 coding sequences, was also independently assembled using reference-guided assembly and V. parahaemolyticus VIPARAQ4037 residues 1748 to 1830 as a reference and was consistently observed in isolates of both the tdh-negative and trh-positive pathotype and the tdh- and trh-positive pathotype and in all instances included trh1 and the T3SS2β genes (Fig. 1B). The association between the trh PAI and T3SS2β with the tdh-negative and trh-positive pathotype and the tdh- and trh-positive pathotype agrees with the findings of previous studies (12, 27). Three clinical tdh- and trh-negative isolates, isolates 04-2548, 09-5357, and 10-4239, did not contain the tdh PAI, the trh PAI, or a T3SS2; however, an almost complete trh PAI including a partial trh gene, urease gene cluster, and the T3SS2β genes was identified in isolate 10-4238, which was categorized as tdh and trh negative by PCR analysis by Banerjee et al. (26) (Fig. 1C). In addition to the TRH hemolysins, the T3SSβ present in this PAI also contains effectors thought to be involved in enterotoxicity and cytotoxicity (5). T3SS2β appears to be a recent acquisition by V. parahaemolyticus and is sometimes found in non-O1, non-O139, and CTX V. cholerae strains (12). The presence of the trh PAI in 10-4238 made us question whether this strain should indeed be classified as a clinical tdh- and trh-negative isolate or if it would be more accurate, on the basis of its WGS, to classify it as a tdh-negative and trh-positive strain. Therefore, this strain was removed from the remaining analysis specific to clinical tdh- and trh-negative isolates.
General genomic features of V. parahaemolyticus clinical isolates.To depict the genetic diversity of pathogenic V. parahaemolyticus strains, 38 clinical isolates (Table 1) representing each of the four previously described genotypes (tdh and trh positive, tdh negative and trh positive, tdh positive and trh negative, and tdh and trh negative) (40, 41) were extensively compared. The pangenome of clinical V. parahaemolyticus isolates was calculated and consisted of 8,399 protein-coding genes (Fig. 2A). To assess the accuracy of computing of a pangenome using draft genomes, three clinical isolates with closed genomes (see Table S1 in the supplemental material) were added to our data set and the pangenome size was recalculated. The addition of three closed genomes caused the size of the pangenome of the clinical isolates to increase to 8,609 genes. This increase would also have been expected after the addition of three draft genomes, and therefore, we concluded that our calculations based on the draft genomes were accurate.
General genomic features of clinical V. parahaemolyticus isolates. (A) The pangenome of clinical V. parahaemolyticus isolates was constructed using the de novo assembly of 38 clinical isolates and contained 8,399 genes. Pangenomes were also assembled for each pathotype and compared. (B) A phylogenetic tree, constructed from concatenated core genes, shows the phylogeny of clinical isolates and demonstrates that each of the pathotypes is polyphyletic.
A gene intersection analysis of the accessory genomes was performed. From this analysis, we defined a core genome of orthologous genes that were shared by all clinical V. parahaemolyticus isolates. The core genome contained 3,807 protein-coding genes, which represented between 76 to 81% of each isolate's genome (Fig. 2A). The size of the core genome remained relatively stable for each additional genome added after the first four. When the closed genomes of three additional clinical isolates (see Table S1 in the supplemental material) were added to the data set and the core genome was recalculated, it decreased to a size of 3,803 protein-coding genes. This indicated that the 38 draft genomes sequenced here provided an excellent estimation of the true core genome of clinical V. parahaemolyticus isolates.
Accessory genes that were unique to each pathotype were also identified, and 654, 520, 918, and 1,097 genes were specific to the tdh- and trh-positive, tdh-negative and trh-positive, tdh-positive and trh-negative, and tdh- and trh-negative pathotypes of clinical isolates, respectively (Fig. 2A). Strain-specific unique genes were also identified, and the sizes of the unique genome varied between strains. For example, seven strains (09-3216, 10-4303, 10-7197, 09-3218, 09-4434, 09-4660, and 09-4664) possessed no unique genes, while one strain, 04-2548 (tdh and trh negative), had 405 unique genes. Clinical isolates of the tdh- and trh-negative pathotype consistently had large unique genomes.
To determine if clinical tdh- and trh-negative isolates are monophyletic and possibly the result of a single loss of a pathogenicity island, an unrooted phylogenetic tree was constructed from the concatenated amino acid sequences of the 3,899 core genes (Fig. 2B). Use of the core genes provides a high-resolution view of phylogeny. This revealed that isolates of each of the pathotypes, tdh- and trh-negative isolates, are polyphyletic. This may indicate a high degree of mobility of pathogenicity elements between V. parahaemolyticus isolates, rather than a single PAI deletion event.
The distribution of clusters of orthologous groups (COGs) of proteins was determined for the core genome, the accessory genome of each pathotype, and the unique genome of each of the tdh- and trh-negative isolates, to determine if there were differences in the proportion of the genome attributable to particular cellular processes in clinical isolates (see Table S2 in the supplemental material). Almost 20% of the core genome was classified as having an unknown function; this proportion dropped to less than 10% in each of the pathotype-specific unique genomes (Fig. 3A and B). The tdh- and trh-negative isolates had a large proportion of genes involved in cell motility relative to the proportion for the other pathotypes (Fig. 3B). Individual tdh- and trh-negative isolates also had large and functionally consistent unique genomes (Fig. 3C). Within this functional category of cell motility, each of the three clinical tdh- and trh-negative isolates contained a methyl-accepting chemotaxis protein not observed in the other pathotypes or in environmental isolates. Cell motility is generally considered to be a factor associated with the ability of V. parahaemolyticus to survive in vivo (5); therefore, the finding of an increase in the number of genes involved in motility in clinical isolates is logical.
COG profiles of the core (A), accessory (B), and unique (C) genomes of clinical tdh- and trh-negative isolates. The numbers at the top of each column denote the number of genes in the unique genome of the corresponding strain listed at the bottom of the column. Abbreviations: neg, negative; pos, positive.
General genomic features of V. parahaemolyticus environmental isolates.To provide a basis for comparison of the clinical tdh- and trh-negative isolates, we sequenced the genomes of 5 environmental tdh- and trh-negative isolates and collected the genomes of an additional 11 from the NCBI database (Table 1). To demonstrate the wide diversity and the relationships of the strains used in this study, a rooted phylogenetic tree was constructed from the concatenated nucleotide sequences of seven housekeeping genes (recA, gyrB, dnaE, dtdS, pntA, pyrC, and tnaA) traditionally used in V. parahaemolyticus MLST analysis (Fig. 4A). Clinical and environmental isolates shared several common lineages, again demonstrating the dynamic nature of virulence factors in this species.
General genomic features of environmental V. parahaemolyticus isolates. (A) A phylogenetic tree, constructed from concatenated MLST sequences, demonstrates the diversity of the strains used in this study, as well as the relationships between each clinical and environmental strain included in this study. (B) A rarefaction curve of the genetic diversity of clinical and environmental V. parahaemolyticus strains was created. Environmental isolates have a much greater genetic diversity than clinical isolates. In addition, these curves demonstrate that the pangenome of V. parahaemolyticus is open.
The sizes of the core genome and pangenome were calculated for the 16 environmental isolates. The core genome of the environmental tdh- and trh-negative isolates was composed of 2,773 protein-coding genes, and though the pangenome was constructed from fewer genomes, it was much larger than that of the clinical isolates at 11,669 protein-coding genes (Fig. 4B). For both the clinical and environmental isolates, the pangenome size increased, in terms of protein-coding gene number, after the addition of each genome, which indicates an open pangenome (42), and this is in agreement with the findings of earlier studies (43) (Fig. 4B). The increased size of the pangenome in environmental isolates is likely due to their ability to survive in diverse niches in situ. For example, environmental isolates were collected from colonized shrimp, muscle, oyster, and clams as well as from marine water and sediments, while clinical isolates have been preselected by use of a very narrow criterion, which is the ability to colonize humans and cause illness. This likely leads to the low level of genetic diversity observed in clinical isolates.
Whole-genome comparison of clinical and environmental tdh- and trh-negative isolates.The genomes of the 3 clinical tdh- and trh-negative isolates were compared with the genomes of 16 environmental tdh- and trh-negative isolates, resulting in the identification of 862 protein-coding genes unique to clinical isolates (Fig. 5; see also Table S3 in the supplemental material). From these genes, 529 were annotated as hypothetical proteins. A large portion of both hypothetical proteins and annotated proteins had high sequence similarity with other genes from other enteric bacteria, such as V. cholerae, Listeria monocytogenes, Campylobacter jejuni, Salmonella enterica, Escherichia coli, and Enterobacter. While this finding is not conclusive, it may indicate that the acquisition of genes from other enteric bacteria may contribute to the ability of clinical tdh- and trh-negative isolates to colonize humans during a gastrointestinal illness.
Genomic comparison of environmental and clinical tdh- and trh-negative isolates. On the basis of a comparison of protein-coding genes, 862 protein-coding genes that are unique to clinical tdh- and trh-negative isolates were identified. neg, negative.
Pathogenicity islands.By examining the complete genome of V. parahaemolyticus RIMD2210633, Hurley et al. (2006) identified seven genomic islands which occur in pathogenic V. parahaemolyticus isolates (38). VPAI-1, VPAI-4, VPAI-5, and VPAI-6 appeared to represent DNA acquired by pandemic V. parahaemolyticus isolates (38). The presence of these islands was variable in our clinical and environmental isolates (Tables 2 and 3). VPAI-1 was present in most clinical tdh-positive and trh-negative isolates but was largely absent from isolates of the other pathotypes. The entire VPAI-2 was present in most tdh-positive and trh-negative isolates and tdh-negative and trh-positive isolates but was only partially present in tdh- and trh-positive and tdh- and trh-negative strains. VPAI-6 was consistently present in each of the clinical isolates. Strains 07-2965 and 08-0278 did not typically have the same PAI profile as the other tdh-positive and trh-negative strains, and this was also reflected in the phylogenetic tree (Fig. 2B). The 10-4238 tdh- and trh-negative strain carried VPAIs more similar to those carried by the tdh-negative and trh-positive isolates, agreeing with findings presented earlier in this paper that this strain is likely a tdh-negative and trh-positive strain and was misidentified by PCR analysis. The environmental tdh- and trh-negative isolates carried pathogenicity island profiles similar to those of the clinical isolates (Table 3), this observation raises several questions about the true roles of these islands and whether their inclusion in pandemic strains is associated more with fitness in the environment than with pathogenicity.
Distribution of pathogenicity islands in clinical isolates
Distribution of pathogenicity islands in environmental isolates
IslandViewer3 was used to search the closed genome of CDC_K4557, which is also a clinical V. parahaemolyticus tdh- and trh-negative isolate, for genomic islands that are common to clinical isolates but that are not found in environmental tdh- and trh-negative isolates (37). We reasoned that if an island was present in at least some of the clinical tdh- and trh-negative isolates but absent from all of the environmental tdh- and trh-negative isolates, it would be a reasonable candidate for evaluation as a novel pathogenicity island. IslandViewer3 was used to identify 29 genomic islands on chromosome 1 and 8 genomic islands on chromosome 2 (Fig. 6). These genomic islands were assessed for their presence or absence across our 43 genomes (see Table S4 in the supplemental material). We found that islands 1, 3, 33, 35, and 36 (as denoted in Fig. 6) were present in almost every strain of V. parahaemolyticus. An island that was present in clinical tdh- and trh-negative strains but not in environmental tdh- and trh-negative strains was not identified.
Genomic islands identified in CDC_K4557. A closed genome is required to identify novel genomic islands. To determine if a novel pathogenicity island is responsible for virulence in the tdh- and trh-negative pathotype, the closed genome of V. parahaemolyticus CDC_K4557, another clinical tdh- and trh-negative isolate, was searched for genomic islands using IslandViewer3, and the images shown above were modified from the IslandViewer3 (37) website (http://www.pathogenomics.sfu.ca/islandviewer/accession/NC_021848.1/ [top panel] and http://www.pathogenomics.sfu.ca/islandviewer/accession/NC_021822.1/ [bottom panel]). BLASTn was used to search for each of the islands in the clinical and environmental isolates from Table 1. IslandViewer3 uses multiple algorithms (IslandPick, SIGI-HMM, and IslandPath-DIMOB) to predict the presence of genomic islands. Red, islands identified by an integrative algorithm incorporating multiple mechanisms of island prediction; blue, islands predicted by the IslandPath algorithm, yellow, islands predicted by the SIGI-HMM algorithm; green, islands predicted by the IslandPick algorithm. M, millions; k, thousands.
Type VI secretion systems.Two T6SSs have previously been found in the pangenome of V. parahaemolyticus. T6SS1 (VP1387 to VP1414) is found on chromosome 1 and is commonly associated with clinical isolates, while T6SS2 (VPA1025 to VPA1046) is found on chromosome 2 and has been found in all tested strains (16, 20). In the current study, six clinical isolates (10-4238, 08-0278, 09-5357, 07-2965, 09-3217, and 10-4255) and eight environmental isolates (ISF-54-12, S171, VPCR-2009, VPTS-2009, SG176, 22720, J-C2-34, and VH3) did not have a T6SS1. Of the isolates that had a T6SS1, variance in this gene cluster was observed between isolates, although variance occurred only in the VP1388, VP1389, and VP1390 genes (Fig. 7). There were five different alleles of these genes, and a phylogenetic tree is shown to demonstrate which isolates contained which alleles (Fig. 7). The 10-4239 isolate had unique alleles for each of these genes which were not observed in any of the other isolates. The N terminus of VP1388 was conserved in all isolates that had a T6SS1, while variance was observed in the C terminus. There is biological significance underlying the variation observed in VP1388 to VP1390. VP1388 has previously been identified to be an antibacterial effector, and VP1389 is its associated immunity protein (21). Changes in the effector must be accompanied by changes in the immunity protein to maintain self-protection. The putative functionality of the third larger gene downstream (VP1390) is still unknown, but on the basis of its association with VP1388 and VP1390, it may have a role in antimicrobial activity. The finding of a novel putative effector in a V. parahaemolyticus clinical tdh- and trh-negative isolate indicates that this protein should be further investigated for roles during infection.
T6SS1. The gene cluster homologous to T6SS1 is present in 32 of our clinical isolates; however, variance in the VP1388 to VP1390 genes was observed between isolates.
Conclusion.The ability of tdh- and trh-negative strains to cause clinical illness is still controversial, and several theories have been proposed to explain why tdh- and trh-negative strains are sometimes isolated from clinical cases, including coinfection with pathogenic V. parahaemolyticus strains, the loss of virulence genes during infection, the presence of novel and uncharacterized virulence factors, or the fact that they play a role in a multistrain infection. However, in this investigation we have identified 862 genes that are present in clinical tdh- and trh-negative isolates but that are not present in environmental isolates. Several of these genes are highly homologous to genes from other enteric bacteria, indicating that horizontal gene transfer may play an important role in the ability of tdh- and trh-negative isolates to survive in the human gastrointestinal tract. In addition, tdh- and trh-negative isolate 10-4239 contains a unique T6SS1 effector/immunity gene combination that should be investigated further.
ACKNOWLEDGMENTS
We are very grateful to Robyn Kenwell for technical assistance with the sequencing of the genomes analyzed in this paper. We also thank John Austin and Franco Pagotto of the BMH Research Division of Health Canada for reviewing the manuscript and offering helpful comments.
This work was funded (A-base) by Health Canada in support of Canada's Food Safety Programs. J.R. is supported by the Visiting Fellow in a Government Laboratory program.
J.R., N.P., C.C.L., and S.K.B. conceived of and designed the experiments. J.R. and C.C.L. performed all wet laboratory experimental, molecular, and next-generation sequencing work. J.R. and N.P. conducted the in silico analysis, including the use of novel tools contributed by A.W.P.; J.R. wrote the paper with input from all other authors.
FOOTNOTES
- Received 22 October 2015.
- Accepted 25 November 2015.
- Accepted manuscript posted online 4 December 2015.
Supplemental material for this article may be found at http://dx.doi.org/10.1128/AEM.03465-15.
- Copyright © 2016, American Society for Microbiology. All Rights Reserved.
