|Home | About | Journals | Submit | Contact Us | Français|
Thirty-two genome sequences of various Vibrionaceae members are compared, with emphasis on what makes V. cholerae unique. As few as 1,000 gene families are conserved across all the Vibrionaceae genomes analysed; this fraction roughly doubles for gene families conserved within the species V. cholerae. Of these, approximately 200 gene families that cluster on various locations of the genome are not found in other sequenced Vibrionaceae; these are possibly unique to the V. cholerae species. By comparing gene family content of the analysed genomes, the relatedness to a particular species is identified for two unspeciated genomes. Conversely, two genomes presumably belonging to the same species have suspiciously dissimilar gene family content. We are able to identify a number of genes that are conserved in, and unique to, V. cholerae. Some of these genes may be crucial to the niche adaptation of this species.
The species concept for bacteria has long been under siege from several angles, and now with thousands of bacterial genomes being sequenced, the disputes have intensified . One frequently used definition of a bacterial species is “a category that circumscribes a (preferably) genomically coherent group of individual isolates/strains sharing a high degree of similarity in (many) independent features, comparatively tested under highly standardized conditions” . Such independent features are usually phenotypes that can easily be tested. For a new species to be defined, amongst other criteria, inter-species DNA–DNA hybridisation has to be below 70%, although this rule is not without its limitations . In the late 1970s and 1980s, the 16S rRNA gene sequence was introduced as a molecular clock that could be used to infer phylogenetic relationships . Ideally, isolates belonging to the same species have identical or nearly identical 16S rRNA genes, and these differ from isolates belonging to different species [32, 44]. In practice, this is not always the case. Examples exist of different species sharing identical rRNA genes (for instance, E. coli and Shigella  that are even placed in different genera); in addition, isolates of one species can have different rRNA genes beyond the 97% that is considered to demarcate species . Lateral transfer of genetic material (to which ribosomal genes are believed to be resistant) destroys the phylogenetic relationship, so that phylogenies based on alternative housekeeping genes can differ from a 16S rRNA tree and frequently are not even in accordance to each other. Such observations question the validity of a phylogenetic tree as the most suitable model for bacterial ancestry, when multiple genetic transfers would produce a network-like evolutionary structure . On the other hand, it is observed that lateral gene transfer is most frequent between genetically related members sharing a similar base content and occupying the same ecological niche . Nevertheless, a core of genes can be recognised that produce coherent phylogenetic trees, though these may not represent the species’ complete evolutionary history as they comprise only a minor fraction of the genetic content of the organism .
Whether a tree or a network is more accurate to describe phylogeny, in either case bacterial species may be considered as a cloud of isolates having a higher level of genetic similarity to each other than to organisms belonging to a different species. When such clouds have fuzzy and overlapping borders, the species concept falls apart but that will only apply to certain cases . Since 16S rRNA genes are not informative on the level of diversity within a species, the 'density' of a cloud of isolates making up a species cannot be determined by this gene. Those genes shared by all isolates belonging to one species comprise the core genome of that species , and the degree of diversity in the remaining non-core genes determines the density of the species cloud.
We hypothesised that certain genes can be recognised as specific to a particular species, to be conserved in that species but not present in related species. We tested our hypothesis with complete genome sequences of the bacterial family Vibrionaceae, which belong to the γ-Proteobacteria and comprises eight genera. Most available genome sequences belong to the genus Vibrio. This genus contains 51 recognised species [10, 46] which are mainly found in marine environments, frequently living in association with marine organisms such as corals, fish, squid or zooplankton. Most of them are symbionts and only a few are human pathogens, notably particular serotypes of V. cholerae producing cholera, Vibrio parahaemolyticus (causing gastroenteritis) and Vi vulnificus (causing wound infections) . Other Vibrionaceae, including V. vulnificus, Aliivibrio salmonicida and V. harveyi, are fish or shellfish pathogens and have major economic impact. Photobacterium profundum, representing another genus within the Vibrionaceae, was also included.
The gene content of 32 available sequenced Vibrionaceae genomes was compared and the results were analysed in various ways. The data allowed us to identify possible V. cholerae-specific genes, since this species was represented by 18 genomes that was a sufficient number to test conservation both within the species and across species. We found that a two-component signal transduction pathway is uniquely conserved in V. cholerae but is not found outside this species. Our findings further indicated that possibly a relatively small set of genes could confer niche specialisation allowing V. cholerae to be adopted to a unique environment, so that over time V. cholerae have become a distinct species.
Publicly available genome sequences of Vibrionaceae were selected that were provided in less than 300 contigs and in which full-length 16S rRNA sequence could be found using the rRNA gene finder RNAmmer . The 32 genome sequences included are shown in Table 1.
The gene annotations as provided in GenBank were used, except for those genomes marked “Easygene” in Table 1 where protein annotation was not available in the RefSeq file at the time of analysis, and we used EasyGene  to identify the genes. As a control, an available GenBank annotation was compared to a generated Easygene annotation to confirm that the number of identified genes was comparable.
RNAmmer  was used to identify 16S rRNA sequences within the 32 genomes. Sequences were considered reliable if they were between 1,400 and 1,700 nucleotides long and had an RNAmmer score above 1,800. In cases where the program found multiple and variable 16S sequences within a genome, one of these (with satisfactory RNAmmer scores) was arbitrarily chosen. The sequences were aligned using PRANK [23, 24], and the program MEGA4 was used to elucidate a phylogenetic tree . Within MEGA4, the tree was created using the Neighbor-Joining method with the uniform rate Jukes–Cantor distance measure and the complete-delete option. Five hundred resamplings were done to find the bootstrap values.
Clustering based on shared gene families from the Vibrio pan-genome was constructed, based on BLASTP similarity using default settings. A BLASTP hit was considered significant if the alignment produced at least 50% identity for at least 50% of the length of the longest gene (either query or subject). Using this criterion, each pair of genes producing a significant reciprocal best hit was scored as belonging to the same gene family. A genome matrix was constructed, containing one row for each genome and one column for each gene family. Cell (i, j) in this matrix is 1 if genome i has a member in gene family j, 0 otherwise. A hierarchical clustering, with average linkage based on the Manhattan distance between genomes was then performed. Two trees were made, one with more weight given to gene families present in most (90%, or between 27 and 30) Vibrio genomes (“stabilome”), and the other with more weight given to gene families present in only a few (two, three, or four) genomes (“mobilome”). Thus, the original Boolean matrix is now scaled differently, depending on the number of genomes in each gene family . For both trees, singletons (families which are only found in one genome) have been excluded.
The results of the BLAST analysis were also used to construct a pan- and core genome plot as follows. Based on clusterings from the pan-genome family tree, an ordered set of genomes was constructed with V. cholerae genomes at the start. For the first chosen genome, all BLAST hits found in the second genome were recorded and the accumulative number of gene families (as defined above) now recognised in total was plotted for the pan-genome. The number of gene families with at least one representative gene in both genomes was plotted for the core genome. A running total is plotted for the pan-genome which increases as more genomes are added, whilst the core genome representing conserved gene families slowly decreases with the addition of more genomes.
The predicted genes of every genome (annotated or found by Easygene) were translated and every gene was compared, by BLASTP against every other genome and its own genome. In the latter case, the hit to self was ignored. The 50/50 rule for BLAST hits as described above was used. If these requirements were met, genes were combined in a gene family. The BLAST results were visualised in a BLAST matrix , which summarises the results of genomic pairwise comparisons and reports, both as percentage and as absolute numbers, the number of reciprocal BLAST hits as a fraction of the total number of gene families found in the two genomes. For easier visual inspection, the cells in the matrix are coloured darker as the fraction of similarity increases. Hits identified within a genome are differently coloured.
BLAST results were also visualised in a BLAST atlas, this time visualising, for all genes in the reference genome V cholerae N16961, their best hit in all other genomes, again with a threshold of 50% identity over at least 50% of the length of the query protein. The atlas displays the hits as they are located in the reference strain . The BLAST scores obtained for each queried gene is plotted, so that conserved and variable regions are located with respect to the reference genome. Note that genes absent in the reference genome are not shown in the lanes of the query genomes.
A phylogenetic tree based on the 16S rRNA gene extracted from the 32 analysed Vibrionaceae genomes is shown in Fig. 1. The 18 V. cholerae genomes build a tight subcluster, quite distanced from the other species. Above this in the figure, another subcluster comprising eight genomes representing at least six species is recognised, and within this cluster the two V. parahaemolyticus genes are not found on the same branch. A third cluster, a bit further removed, includes Aliivibrio fischeri and A. almonidica as well as V. splendidus and Vibrio species MED 222; the gene of Photobacterium profundum is the most distant.
Starting with a database containing the total set of all Vibrio gene families, a profile of matching gene families was constructed for each individual genome. This was stored as a matrix, containing a column for each gene families, and a row for each genome. The rows contain a 0 or 1 representing the presence or absence of the gene family. This matrix was weighted to emphasise either the genes found in most genomes (the “stabilome”) or in only a few genomes (the “mobilome”); from these weighted matrices, clustering of gene families yielded the resulting trees shown in Fig. 2. Shorter distances represent genomes with many gene families in common, and larger distances reflect genomes with fewer gene families in common. As expected, in both trees, genomes from the same species cluster together, whereby the depth of resolution within a species is considerably better than can be seen in the 16S rRNA tree in Fig. 1. Similarity between the unspeciated Vibrio isolate MED222 and V. splendidus is suggested by their close clustering; this is a connection also suggested by others . Note that the unspeciated Vibrio isolate Ex25 and V. parahaemolyticus 2210633 cluster together in the mobilome tree, but are more distant in the stabilome. This implies that the genes shared between these two genomes are less common genes within the Vibrio genomes examined here. As already indicated by the 16S rRNA tree, the two V. parahaemolyticus isolates are quite dissimilar, and appear on separate branches. The Aliivibrio cluster is placed within Vibrio genomes in both the stabilome and the mobilome, as was the case for their 16S rRNA gene. P. profundum is not such an outlier as in the 16S rRNA tree, and in the stabilome. It is even positioned close to the Aliivibrio genomes. Zooming in at the genomes of V. cholerae, a division into two subclusters can be seen; these clusters correspond to environmental vs. clinical isolates (with the exception of V52 in the stabilome).
BLAST results were analysed to construct a pan-genome, which is a hypothetical collection of all the gene families that are found in the investigated genomes . The core genome was constructed from all gene families that were represented at least once in every genome. Thus, the gene families conserved in all genomes represent their core genome; adding the remaining gene families produces the pan-genome. The resulting pan- and core genome plot is shown in Fig. 3. The genomes start with the documented clinical isolates of V. cholerae and then follow the order suggested by the pan-genome family clustering (Fig. 2), although genomes from the same species were kept together (the two V. parahaemolyticus genomes were split in the trees). As more genomes are added in the plot, the number of gene families in the pan-genome (blue line) increases, and the number of conserved gene families (red line) in the core genome decreases, albeit at a lower rate. This is because every genome can add many novel (and frequently different) genes to the pan-genome but only decreases the core genome with a few genes that are absent in that particular strain but that were conserved in the previously analysed genomes. The pan-genome curve increases with a relative steep slope when a novel species is added, as is obvious when a V. parahaemolyticus genome is added after the last V. cholerae. A stable plateau can be seen for the pan-genome of V. cholerae around 6,500 genes. Nevertheless, a small increase occurs when adding V. cholerae 11587; this is caused by the difference between the two subclusters of V. cholerae seen in Fig. 2. V. cholerae strain 2740-80 behaves atypical in all the figures shown; although documented as an environmental isolate, it appears closer to the clinical isolates, in terms of overall genomic properties.
When the first genome of A. fischeri is added, which is not a member of the Vibrio genus, it does not add significantly more novel genes to the pan-genome than Vibrio genomes did. This contrasts with P. profundum which produces a sharp increase in the pan-genome, as does, interestingly, V. shilonii. Note that there are approximately 20,200 total gene families within the 32 sequenced Vibrionaceae genomes, whereas the core genome decreases to approximately 1,000 gene families.
A BLAST matrix provides a visual overview of reciprocal pairwise whole-genome comparisons, as shown in Fig. 4. The stronger a matrix cell is coloured, the more similarity was detected between the gene content of two genomes. As can be seen in the lower right triangle, all V. cholerae genomes are highly similar, with similarity ranging between 64% and 93% for any given pair of genomes. No statistical difference was observed when comparing clinical isolates to environmental isolates. The two A. fischeri and the two V. vulnificus genomes also share a high degree of identity within their species (75% and 67%, respectively), visible at the bottom of the matrix. In contrast, the two V. parahaemolyticus genomes only share 35% identity, which is not higher than the similarity detected between genomes of different species. With 72% similarity, isolate MED222 most closely matches V. splendidus and with 65% isolate EX25 again shares most similarity with V. parahaemolyticus 2210633.
A BLAST atlas was constructed using V. cholerae N16961 (O1, El Tor) as the reference genome, shown in Fig. 5. The best blast hits identified in the query genomes are plotted in the lanes around the reference genome, with different colours for different species. In general, chromosome 1 is more strongly conserved than chromosome 2. A large part of chromosome 2 of N16961 displays very little conservation in the other genomes; this area represents a super integron  that contains the V. cholerae-specific repeat (VCR) sequences, as well as a high number of gene cassettes. The repeat sequences are visible as black boxes in the repeat lane of the reference genome (second inner lane). Although all V. cholerae genomes contain a superintegron, its genes are very diverse between isolates  which explains the lack of blast hits in this region.
Several regions of the atlas have been highlighted. Gaps B, C, D and F on chromosome 1 (indicated in green) contain genes that are conserved in the represented genomes of V. cholerae but not in the other Vibrionaceae. The gaps marked A, E and G indicate regions that are specific to the toxigenic, clinical isolates only. Annotated, V. cholerae-specific genes present in all these regions are listed in Table 2 (hypothetical genes are excluded). Genes specific for toxinogenic V. cholerae identified in gap A include, amongst others, biosynthesis genes for the toxin co-regulated pilus (which is required for transmission of the prophage CTXΦ carrying the enterotoxin genes), as well as genes encoding citrate lyase. Note that the genes in gap A are also found in the environmental isolate V. cholerae 2740-80.
Gap B contains a number of outer membrane protein genes involved in sugar modification that are found in all V. cholerae genomes. Genes from gap C encoding a histidine kinase two-component signal transduction regulatory system are also conserved within the species, as genes in gaps D and F, involved in chemotaxis and possible multidrug resistance.
Gap E, containing genes conserved in toxigenic strains only, holds the prophage CTXΦ that contains the genes encoding cholera enterotoxin subunits A and B; this enterotoxin is responsible for the excessive, watery diarrhoea typical for cholera. Upon binding to target cell GM1 gangliosides, enterotoxin enters the cell and stimulates adenylate cyclase by ADP ribosylation. The resultant increased cyclic AMP levels induce excessive electrolyte movement and sodium plus water secretion . Strain M66-2 is believed to be a precursor of the seventh pandemic V. cholerae that lacks the prophage CTXΦ and the enterotoxin genes . Gap E bears the RTX toxin operon, which encodes a pore-forming cytotoxin . An RTX toxin is also present in environmental isolate 2740-80 and in V. vulnificus.
Gap G on chromosome 2 consists of a set of five genes, all in the same orientation, in a putative operon, flanked by genes on the complimentary strand. This appears to be a remnant of a mobile element, as these genes are flanked by a transposase gene on the 3′ end, and there is a small global repeat on the 5′ end. Only the first two of the five genes have an assigned function, with the first gene being a GMP reductase, and the second a putative DNA methyltransferase. The remaining three genes are hypothetical, but their strikingly strong conservation in all pathogenic strains and complete absence of homologues in the other Vibrio genomes strongly point towards a potential biological significance.
The recent availability of many Vibrionaceae genomes, including a substantial number of V. cholerae genomes, allows the possibility to take a closer look at the similarities and differences of species within the genus Vibrio. This can examine, on a genome scale, what distinguishes V. cholerae from the other Vibrio species. Since not all V. cholerae isolates are pathogenic, the presence of the prophage-bearing cholera enterotoxin, the main virulence factor for cholera, is not a suitable marker for this species. We attempted to identify a set of V. cholerae-specific genes, and also explored the internal diversity within the V. cholerae genomes that have been sequenced to date.
On a phylogenetic tree based on the 16S ribosomal RNA gene, those isolates that do not belong to the genus Vibrio were positioned as outliers, as expected. This tree further indicated the closest resembling 16S rRNA sequence for the two sequenced Vibrio strains that are currently not assigned to a species. It was observed that the two sequenced V. parahaemolyticus strains were not placed together. The complete gene content of each genome was next compared by BLAST and the results were pooled into gene families which were subjected to cluster analysis. This provided evidence that the 18 V. cholerae genomes fall into two subclusters, one mainly containing clinical isolates and the other environmental isolates.
The gene family clustering, subsequent pan-genome analysis and the pairwise BLAST results, as summarised in the BLAST matrix, all supported the relatedness of Vibrio species Ex25 to V. parahaemolyticus 2210633 but not to V. parahaemolyticus 16. This latter genome was quite different from V. parahaemolyticus 2210633 in all analyses. Although it is possible that the species V. parahaemolyticus is far more genetically diverse than V. cholerae, A. fischeri or V. vulnificus, an alternative explanation is that one of the sequenced isolates is perhaps incorrectly named as V. parahaemolyticus. The similarity between Vibrio species MED222 and V. splendidus based on gene families is in agreement with their related 16S rRNA genes and published data . However, in contrast to what the ribosomal gene suggests, our whole-genome comparison indicates that the three Aliivibrio genomes (A. salmonicida and two A. fischeri) are not so different from Vibrio after all. Their recent placement in the genus Aliivibrio, a decision based on five genes (the 16S rRNA gene and four housekeeping genes) and phenotypical characteristics , appears not to be reflective of the whole genome picture presented here.
The BLAST results were graphically summarised in a BLAST atlas, which visualised V. cholerae-specific gene clusters. These coded for polysaccharide biosynthesis enzymes, response regulators and chemotaxis proteins, amongst others. In addition, a V. cholerae-specific, histidine kinase two-component signal transduction regulatory system was identified. The two-component signal transduction pathway is a powerful regulating system for bacteria to adapt to a particular ecological niche. There is a precedent for this claim, as the introduction of a single regulatory protein in Vibrio fischeri strain MJ11 has been shown to specifically enable colonization of the squid Euprymna scolopes .
As expected, the main differences observed between V. cholerae clinical isolates and the environmental strains are due to genes related to virulence. Two exceptions are the presence of a number of virulence genes in the environmental strain V. cholerae 2740-80 and the absence of enterotoxin genes in clinical isolate M66-2. It has already been suggested that M66-2 might be a predecessor of pandemic, enterotoxic V. cholerae . From sequence comparison of four housekeeping genes, it was concluded that V. cholerae 2740-80 is intermediary between toxigenic and non-toxigenic isolates . This view is confirmed by the data presented here, although we propose to consider the possibility that the isolate arose from a pandemic clone that has lost the CTXΦ prophage, rather than being a precursor of a pathogen.
In conclusion, several different methods of genome comparisons have yielded a picture of V. cholerae genomes as forming a distinct cluster, compared to related species, and a relatively small number of genes might be responsible for environmental niche adaptation and hence for generation of this distinct species. Likely candidates include multiple two-component signal transduction regulatory proteins as well as chemotaxis proteins.
We would like to thank Tim Binnewies for early work on this project, and also to the Danish Research Councils and the DTU Globalization funds for financial support.
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.