Carbohydrate metabolic systems present on genomic islands are lost and gained in Vibrio parahaemolyticus

Abstract

Background: Utilizing unique carbohydrates or utilizing them more efficiently help bacteria expand and colonize new niches. Horizontal gene transfer (HGT) of catabolic systems is a powerful mechanism by which bacteria can acquire new phenotypic traits that can increase survival and fitness in different niches. In this work, we examined carbon catabolism diversity among Vibrio parahaemolyticus, a marine species that is also an important human and fish pathogen.

Results: Phenotypic differences in carbon utilization between Vibrio parahaemolyticus strains lead us to examine genotypic differences in this species and the family Vibrionaceae in general. Bioinformatics analysis showed that the ability to utilize D-galactose was present in all V. parahaemolyticus but at least two distinct transporters were present; a major facilitator superfamily (MFS) transporter and a sodium/galactose transporter (SGLT). Growth and genetic analyses demonstrated that SGLT was a more efficient transporter of D-galactose and was the predominant type among strains. Phylogenetic analysis showed that D-galactose gene galM was acquired multiples times within the family Vibrionaceae and was transferred between distantly related species. The ability to utilize D-gluconate was universal within the species. Deletion of eda (VP0065), which encodes aldolase, a key enzyme in the Entner-Doudoroff (ED) pathway, reached a similar biomass to wild type when grown on D-gluconate as a sole carbon source. Two additional eda genes were identified, VPA1708 (eda2) associated with a D-glucuronate cluster and VPA0083 (eda3) that clustered with an oligogalacturonide (OGA) metabolism cluster. EDA2 and EDA3 were variably distributed among the species. A metabolic island was identified that contained citrate fermentation, L-rhamnose and OGA metabolism clusters as well as a CRISPR-Cas system. Phylogenetic analysis showed that CitF and RhaA had a limited distribution among V. parahaemolyticus, and RhaA was acquired at least three times. Within V. parahaemolyticus, two different regions contained the gene for L-arabinose catabolism and most strains had the ability to catabolism this sugar.

Conclusion: Our data suggest that horizontal transfer of metabolic systems among Vibrionaceae is an important source of metabolic diversity. This work identified four EDA homologues suggesting that the ED pathway plays a significant role in metabolism. We describe previously uncharacterized metabolism islands that were hotspots for the gain and loss of functional modules likely mediated by transposons.

Keywords: CRISPR-Cas systems; Citrate fermentation; Entner-Doudoroff aldolase (EDA); L-arabinose; L-rhamnose; Metabolism islands; Sodium/galactose transporter SGLT; Tn7-like transposon.

Fig. 1

d-galactose utilization gene cluster. a. Comparative genomic analysis of the d-galactose catabolism and transport region from RIMD2210633, UCM-V493, V. cholerae N16961 and E. coli MG1655. Gray shade, region of nucleotide homology. b. d-galactose utilization pathway showing enzymes and open reading frames (ORFs) designations in RIMD2210633 and UCM-V493. c. RIMD2210633 (open circles) and UCM-V493 (open squares) were grown aerobically at 37 °C for 48 h in M9 + 10 mM d-galactose. d. Growth analysis of RIMD2210633, UCM-V493 and RIMDpSGLT (the sglt gene (ORF VPUCM_0844) was cloned into RIMD2210633) in M9 + 10 mM d-galactose after 24 h. e. Growth analysis of UCM-V493 and Δsglt at 37 °C in M9 + 10 mM d-galactose

Fig. 2

Phylogenetic analysis of a. GalM, b. MFS, and c. SGLT. Homologues of GalM, SGLT and MFS proteins from members of the Vibrionaceae were obtained from the NCBI genome database. The evolutionary history of each protein was inferred using the Neighbor-Joining method [26]. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) are shown next to the branches [27]. Each tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the Dayhoff matrix based method [28] and are in the units of the number of amino acid substitutions per site. The rate variation among sites was modeled with a gamma distribution (shape parameter = 5). All ambiguous positions were removed for each sequence pair (pairwise deletion option). Evolutionary analyses were conducted in MEGA X [29]

Fig. 3

Three keto-deoxy-phosphogluconate aldolases (EDA1, EDA2, and EDA3). a. Gene cluster of d-gluconate, d-glucuronate and oligogalacturonide (OGA) present in V. parahaemolyticus RIMD2210633. b. Predicted catabolic pathway for the utilization of OGA, d-glucuronate and d-gluconate based on the ORFs present in V. parahaemolyticus RIMD2210633

Fig. 4

d-gluoconate metabolism. a. Growth analysis of V. parahaemolyticus RIMD2210633 and Δeda1 mutant strain in M9 supplemented with 10 mM D-glucose and 10 mM D-gluconate. b. Expression analysis of the first gene in the pentose phosphate pathway (VP1708) and the three putative aldolases (VP0065, VPA0083, VPA1708) in wild-type strain in gluconate relative to the expression of these genes in wild-type in glucose. Expression analysis of VP1708 and the two putative aldolases (VPA0083, VPA1708) in the Δeda1 mutant strain in gluconate relative to the expression of these genes in the wild-type in gluconate. *p < 0.05, **p < 0.01, ***p < 0.001

Fig. 5

Phylogenetic analysis of EDA1, EDA2, EDA3 and EDA4. This analysis involved 108 amino acid sequences representing EDA1, EDA2, EDA3 and EDA4 homologues present in the Vibrionaceae. The evolutionary history was inferred using the Neighbor-Joining method [26]. The optimal tree with the sum of branch length = 6.21447650 is shown. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) are shown next to the branches [27]. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the Dayhoff matrix based method [28] and are in the units of the number of amino acid substitutions per site. The rate variation among sites was modeled with a gamma distribution (shape parameter = 5). All ambiguous positions were removed for each sequence pair (pairwise deletion option). There were a total of 213 positions in the final dataset. Evolutionary analyses were conducted in MEGA X [29]

Fig. 6

A 73-kb metabolic island in V. parahaemolyticus. a. Detailed analysis of 73-kb metabolic island identified in chromosome 2 of CDC_K4557, which was missing from RIMD2210633. Arrows indicate ORFs and ORFs with identical color indicate similar function. Gray arrows, genes coding hypothetical and other functional proteins. Black arrow represents a transposase (Tnp). b. Detailed analysis of the citrate fermentation cluster in CDC_K4557 and comparison of the same region in V. cholerae. c. Citrate fermentation in Simmons citrate slant. Conversion of green slant to blue indicates citrate utilization

Fig. 7

Phylogenetic analysis of CitF among the Vibrionaceae. This analysis involved 32 amino acid sequences representing CitF homologues present among the Vibrionaceae. Phylogeny was inferred by using the Maximum Likelihood method and Le_Gascuel_2008 model [40]. The tree with the highest log likelihood (− 7580.58) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model, and then selecting the topology with superior log likelihood value. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 0.4398)). The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. All positions with less than 95% site coverage were eliminated, i.e., fewer than 5% alignment gaps, missing data, and ambiguous bases were allowed at any position (partial deletion option). There were a total of 489 positions in the final dataset. Evolutionary analyses were conducted in MEGA X [29]

Fig. 8

A 135-kb Metabolism island in V. parahaemolyticus. a. Detailed analysis of a 135-kb metabolic island identified in chromosome 2 of FORC_022 absent from RIMD2210633. Arrows indicate ORFs and ORFs with identical color indicate similar function. Gray arrows, genes coding hypothetical and other functional proteins. Black arrows represent transposases. b. Transposases and direct repeats identified flanking the functional modules with the 135-kb island

Fig. 9

Analysis of the 26-kb l-rhamnose utilization island. a. Gene cluster of the l-rhamnose utilization in V. parahaemolyticus AQ3810. b. Genomic locus of the 26-kb l-rhamnose region in AQ3810. c. A Tn7-like transposon and a mini CRISPR-Cas system associated with T3SS-2α and a transposon like region associated with the 26-kb L-rhamnose region. d. Schematic of the L- rhamnose genomic loci in AQ4037, S03-S05, S08-S13 VPA1309 in chromosome 2 that is also the location of T3SS-2α in RIMD2210633. Location of the 26-kb l-rhamnose island in FORC_023 and the empty site at VPA1309. Parenthesis indicates homologous ORFs. e. Phylogenetic analysis of RhaA protein from FORC_022 and AQ3810 among members of the family Vibrionaceae. The evolutionary history was inferred by using the Maximum Likelihood method and Le_Gascuel_2008 model [40]. The tree with the highest log likelihood (− 6340.56) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model, and then selecting the topology with superior log likelihood value. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 0.4168)). The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. This analysis involved 29 amino acid sequences. All positions with less than 95% site coverage were eliminated, i.e., fewer than 5% alignment gaps, missing data, and ambiguous bases were allowed at any position (partial deletion option). There were a total of 418 positions in the final dataset. Evolutionary analyses were conducted in MEGA X [29]