Comparative genomic and transcriptomic analyses reveal distinct response strategies to hypoxia by Vibrio parahaemolyticus isolates of clinical and aquatic animal origins

Abstract

Purpose

Vibrio parahaemolyticus is a leading seafood borne pathogen worldwide. The aim of this study was to decipher the response mechanism of V. parahaemolyticus isolates of clinical and aquatic animal origins to the hypoxic condition, which challenges the bacterial survival in the host and in the environment.

Methods

Growth profiles of V. parahaemolyticus isolates (n = 5) of clinical and aquatic animal origins were examined at different stress conditions (osmolality, acid, temperature, and O₂ concentrations). Draft genomes of the V. parahaemolyticus isolates were determined using the Illumina sequencing technique. Comparative genomic analysis were performed to identify and validate the hypoxic tolerance-related genes.

Results

The V. parahaemolyticus isolates had an oxygen concentration-dependent growth mode, and the 10% O₂ condition strongly inhibited the bacterial growth, when incubated in TSB medium (pH 8.5, 3% NaCl) at 37 °C. Unexpectedly, in marked contrast to the normal 21% O₂ condition, the 10% O₂ treatment for 24 h significantly increased biofilm formation of V. parahaemolyticus isolates (p < 0.05). Draft genome sequences of four V. parahaemolyticus isolates of aquatic animal origins were determined (4.914–5.3530 Mb), which carried mobile genetic elements (n = 12–29). Genome-wide gene expression changes triggered by the hypoxic condition were further examined. Comparative transcriptomic analyses unveiled multiple molecular strategies employed by the bacterium to mitigate the cell damage caused by the hypoxia. Of note, the pathogenic V. parahaemolyticus ATCC17802 down-regulated and/or shut down ten metabolic pathways to reduce cell viability and maintain cell structure under the hypoxic stress.

Conclusions

The results of this study fill prior gaps in the response mechanism of V. parahaemolyticus to the hypoxic condition. Different tolerance to hypoxia contributes to the persistence of pathogenic V. parahaemolyticus in the niches.

Introduction

Vibrio parahaemolyticus is a Gram-negative bacterium that inhabits in marine and estuarine environments worldwide (Li et al. 2019). Consuming raw seafood, incomplete food processing, or cross-contamination during food processing can lead to V. parahaemolyticus infection in humans. The clinical manifestations are gastroenteritis, organ infection, septicaemia, and even death (Li et al. 2019). In the United States, more than 50% of foodborne gastrointestinal-Vibriosis cases are caused by V. parahaemolyticus (Karan et al. 2021), which leads to about 80,000 illnesses each year (https://www.cdc.gov/Vibrio/, accessed on November 1, 2022). In China, a national surveillance of 152,792 patients of all ages with acute diarrhea in 31 provinces was administered in 2009–2018. V. parahaemolyticus was found to be the third common bacterial pathogen and contributed to 10.83% of all positive detection (Wang et al. 2021). The pathogenicity of V. parahaemolyticus is strongly associated with the expression of thermostable direct hemolysin (TDH) and/or TDH-related hemolysin (TRH) (Pazhani et al. 2021).

V. parahaemolyticus is frequently detected in seafood (Vu et al. 2022). The elevated seawater temperatures caused by rising global temperatures from climate change may affect transmission of waterborne pathogenic bacteria (Murray et al. 2020; Tong et al. 2022). For instance, the number of days per year suitable for Vibrio in the Baltic Sea and along the north east coast of the United States reached 107 in 2018, double the early 1980s baseline (Murray et al. 2020). As the sea temperature rises, oxygen (O₂) becomes less soluble via a process known as deoxygenation (Gao et al. 2012). Low oxygen condition refers to oxygen partial pressure of 0–10% (Chen et al. 2007). Compared to strictly aerobic or strictly anaerobic conditions, the hypoxic zone has unique biological properties in the environment (Pohl et al. 2022).

The respiration of oxygen is the main source of energy for biological cells (Bueno et al. 2020). The low oxygen condition impacts on cellular enzymatic activity, such as dehydrogenase and superoxide dismutase (Obbard et al. 1994; Reaney et al. 2005), and also change how quickly organisms metabolize and respire (Oschlies et al. 2018). Studies have indicated that waterborne pathogen Vibrio cholerae can exchange common oxygen respiratory terminal oxidases with high-oxygen-affinity terminal oxidases (H–O-ATOs), using the fermentation pathway, with some electron acceptors (EAs) such as fumarate, and nitrate (Bueno et al. 2018; Heidelberg et al. 2000). It has also been reported that biofilm formation of common pathogens Pseudomonas aeruginosa, Streptococcus aureus, and Escherichia coli induces the production of hypoxia-related enzymes such as H–O-ATOs, nitrate reductases, and anaerobic ribonucleotide reductases (Crespo et al. 2016; Jo et al. 2017; Létoffé et al. 2017). Biofilm, constructed primarily by autogenic extracellular polymeric substances, can protect bacteria from the environmental stress, hinder phagocytosis, and thus confer the capacity for long-term colonization and persistence in the host (Thi et al. 2020). The enteric pathogen V. parahaemolyticus can access, colonize, and proliferate within the human intestine with low oxygen level, causing severe diarrhea disease (Bueno et al. 2020; Sun et al. 2022). Nevertheless, to the best of our knowledge, how V. parahaemolyticus responses to the low oxygen condition has not yet been unveiled so far.

In our previous studies, a number of V. parahaemolyticus isolates of aquatic animal origins were isolated, identified and characterized (Su and Chen 2020; Sun et al. 2014; Xu et al. 2022a; Xu et al. 2022b; Yang et al. 2020; Yu et al. 2022; Zhu et al. 2020). Based on these studies, in the present study, we deciphered for the first time the fitness mechanism of V. parahaemolyticus of aquatic animal and human clinical origins under the hypoxic condition. The major objectives of this study were (1) to investigate growth features of V. parahaemolyticus isolates (n = 5) at different concentrations of oxygen (21–5% O₂); (2) to monitor dynamic process of biofilm formation, and observe morphological cell structure changes of the V. parahaemolyticus isolates at the hypoxic condition (e.g., 10% O₂); and (3) to decipher the molecular mechanism of V. parahaemolyticus response to the hypoxic condition by comparative genomics and transcriptomics analyses. The results of this study provide the first experimental evidence for multiple molecular strategies developed by V. parahaemolyticus to deal with the hypoxia.

Materials and methods

V. parahaemolyticus strains and culture conditions

The V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 strains (Table 1) were individually inoculated (1%, v/v) into 5-mL Tryptone Soya Broth (TSB) (Beijing Land Bridge Technology, Beijing, China), and routinely incubated aerobically with shaking at 180 r/min at 37 °C incubator (Shanghai Zhichu Instrument, Shanghai, China). V. parahaemolyticus ATCC17802 was used as a positive control strain, which was first isolated in 1950, leading to the outbreak of human acute gastroenteritis in Japan (Fujino et al. 1953). The overnight cultures were re-inoculated in fresh TSB, incubated to middle logarithmic growth phase (mid-LGP), and used in the following assays.

Table 1 The phenotype features of the V. parahaemolyticus isolates used in this study

Antibiotic resistance and heavy metal tolerance assays

The V. parahaemolyticus isolates was tested for antibiotic susceptibility using the standard disc diffusion method of Clinical and Laboratory Standards Institute (CLSI, M100-S28, 2018), USA. The Mueller–Hinton (MH) medium and ten antibiotic discs were purchased from OXOID, Basingstoke, UK. Escherichia coli ATCC25922 was used as a quality control strain (Su and Chen 2020).

Tolerance of the V. parahaemolyticus isolates to heavy metals was examined using the standard broth dilution testing (microdilution) (Su and Chen 2020; Xu et al. 2022a; Yu et al. 2022). Eight heavy metals (Analytical Reagents) were purchased from Sinopharm Chemical Reagent Co., Ltd., Shanghai, China. E. coli K12 was used as a quality control strain (Su and Chen 2020).

Growth curve assay of the V. parahaemolyticus isolates at different NaCl concentrations, pH and temperatures

The TSB medium was individually adjusted to different NaCl concentrations (0.5%, 1%, 2%, 3%, 4%, and 5%) and pH values (6.0, 6.5, 7.0, 7.5, 8.0, and 8.5) (Sun et al. 2014). Growth curves of the V. parahaemolyticus isolates under the different NaCl (0.5–5%), and pH (6.0–8.5) conditions were individually measured at 37 °C for 40 h using Multimode Microplate Reader (BioTek Instruments, USA). The absorbance value at 600 nm (OD₆₀₀) was used as a parameter for bacterial biomass (Yao et al. 2020). Growth curves of the V. parahaemolyticus isolates in the TSB (pH8.5, 3% NaCl) were also measured at 37 °C, and 25 °C for 40 h, respectively.

Growth curve assay of the V. parahaemolyticus isolates at different concentrations of O₂

Growth curves of the V. parahaemolyticus isolates at different concentrations of O₂ were measured as described previously (Tian et al. 2020) with minor modifications. Briefly, the V. parahaemolyticus isolates were individually inoculated in the 5-mL TSB (pH8.5, 3% NaCl), which was adjusted and equilibrated to different concentrations of O₂ (21%, 18%, 15%, 10%, and 5% O₂) using the nitrogen-blowing instrument (Wen Dong Chemical, Shanghai, China). The O₂ concentrations were monitored using CY-12C Oxygen Meter (Hangzhou Jiachang Electronic Technology Co., Ltd, Hangzhou, China). Under the different O₂ conditions (21–5% O₂), V. parahaemolyticus isolate was individually inoculated into thirty glass culture tubes (15 × 1.6 cm), which were sealed with sealing film (Beideng Medical Co., Ltd., Jiangsu, China), and incubated with shaking at 180 r/min at 37 °C. Instead of a bioreactor, the simple continuous culture condition allowed to measure OD₆₀₀ values every 4 h for 40 h from three independent tubes using the Multimode Microplate Reader (BioTek Instruments, USA). Meanwhile, O₂ concentrations in air within the tubes were measured correspondingly. Low oxygen condition refers to oxygen partial pressure of 0–10%, which corresponds to a range of dissolved oxygen in water body at room temperature and pressure of 0–4.5 mg/L (Chen et al. 2007).

Biofilm formation assay

Biofilm formation of the V. parahaemolyticus isolates was determined using the crystal violet staining method as described in our previous report (Yang et al. 2020). Briefly, 1 mL/well of each V. parahaemolyticus culture (adjusted OD₆₀₀ = 0.4) was inoculated into 12-well bacterial culture plates (Shanghai Sangon Biological Engineeing Technology and Service Co., Ltd., Shanghai, China). The plates were statically incubated at the 21% O₂, and 10% O₂ conditions, respectively. Biofilms formed at 12 h, 24 h, 36 h, 48 h, and 60 h at 37 °C were individually stained, washed, fixed, and measured as described previously (Yang et al. 2020). The 0.1 M phosphate buffer saline (PBS) (pH 7.2–7.4), and 0.25% crystal violet were purchased from the Sangon (Shanghai, China).

Scanning electron microscope (SEM) assay

Preparation of samples for SEM assay was performed according to the method described previously (Yang et al. 2020). Briefly, the V. parahaemolyticus isolates were individually incubated in the TSB (3% NaCl, pH 8.5) at the 21% O₂, and 10% O₂ conditions at 37 °C for 24 h, respectively. Aliquots of 1.5-mL bacterial culture were centrifuged at 8000 g for 5 min at 4 °C, and the cell pellet was wash with 1 × PBS (pH 7.2–7.4, Sangon, China), fixed with glutaraldehyde (Bioengineering Co., Ltd., Shanghai, China), and observed using the SU5000 SEM (Hitachi, Tokyo, Japan, 5.0 kV, × 30,000).

Genome sequencing, assembly, and annotation

The V. parahaemolyticus isolates were routinely incubated in the TSB (3% NaCl, pH 8.5) at 37 °C to the mid-LGP. The genomic DNA was prepared using TIANamp Bacteria DNA Kit (Tiangen Biochemical Technology Co Ltd., Beijing, China). Three separately produced DNA samples were used for each of the V. parahaemolyticus isolates. Only high quality genomic DNA samples (A₂₆₀/₂₈₀ = 1.8–2.0) were subjected to genome sequencing (Xu et al. 2022a).

Whole-genome sequencing was conducted by Shanghai Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) using the Illumina Hiseq × 10 (Illumina, San Diego, CA, USA) platform (Xu et al. 2022a). The insert size of PE150 (pair-end) sequencing was 400 bp. Low-quality sequence filtering; high-quality sequence assembly; annotation of coding sequences (CDSs), rRNA genes, and tRNA genes; and prediction of Clusters of Orthologous Groups (COG) of proteins were performed using the same software with default parameters as described in our recent report (Xu et al. 2022a).

The virulence factor database (http://www.mgc.ac.cn/VFs) and antibiotic resistance gene database (http://arpcard.Mcmaster.ca) were used to detect virulence- and antibiotic resistance-related genes in the V. parahaemolyticus genomes, respectively. Genomic isolands (GIs), prophages, integrons (INs), insertion sequences (ISs), and clustered regularly interspaced palindromic repeats (CRISPR)-Cas systems were predicted using the same software with default parameters as described in our recent report (Xu et al. 2022a) The MGEs were ordered according to the numbers of the Scaffolds on which they were present (see the corresponding Supplementary Tables).

Phylogenetic tree assay

A total of seventy-eight V. parahaemolyticus isolates were subjected for a phylogenetic tree, of which complete genome sequences of seventy-three V. parahaemolyticus isolates were downloaded from the GenBank database, together with draft genome sequence of V. parahaemolyticus ATCC17802 (Table S1). Amino acid data sets of single-copy orthologs of the V. parahaemolyticus genomes, coupled with the four genomes determined in this study, were analyzed using OrthoFinder (version 2.2.6) software (Emms and Kelly 2019). The phylogenetic tree was constructed using the RAxML (version 8) software (Stamatakis 2014) with 1,000 bootstrap replications and a cut-off threshold of ≥ 50% bootstrap values.

Illumina RNA sequencing and analysis

The V. parahaemolyticus isolates were individually incubated in the TSB (3% NaCl, pH 8.5) under the 10% condition at 37 °C for 24 h. Total RNA was extracted using RNeasy Protect Bacteria Mini Kit, QIAGEN RNeasy Mini Kit. The DNA was removed from the extracted RNA samples using RNase-Free DNase Set (QIAGEN Biotech Co. Ltd., Frankfurt, Germany). The Illumina RNA-sequencing was conducted by Shanghai Majorbio Bio-pharm Technology Co. Ltd. (Shanghai, China) using Illumina HiSeq 2500 platform (Illumina, USA). Three separately prepared RNA samples were used for each of the V. parahaemolyticus isolates. The V. parahaemolyticus isolates grown under the 21% O₂ condition were individually used as controls.

Gene expression was analyzed using the RNA-Seq by Expectation–Maximization (RSEM) software (version1.3.3, http://deweylab.github.io/RSEM/). The criteria of fold changes ≥ 2.0 or ≤ 0.5 and p-values < 0.05 were used to define differentially expressed genes (DEGs). Gene set enrichment analyses (GSEA) were performed if the enrichment test p-values were less than 0.05 (Yang et al. 2020; Yu et al. 2022).

Reverse transcription real time-quantitative PCR (RT-qPCR) assay

Reverse transcription was performed using PrimeScriptTM RT reagent Kit with gDNA Eraser (Perfect Real Time) (QIAGEN Biotech Co. Ltd., Hilden, Germany). The TB Green® Premix Ex Taq™ II (TliRNaseH Plus) (Takara Biomedical Technology Co., Ltd., Beijing, China) was used for relative quantitative PCR with 16 s RNA as the internal reference gene. The RT-qPCR reaction was performed using the Fast Real-Time analyzer (Applied Biosystems, Foster City, California, USA). with the following protocol: 50 °C for 2 min; 95 °C pre-denaturation for 10 min; and 95 °C denaturation for 15 s, 60 °C annealing for 1 min, for 40 cycles. The relative expression levels of target gene and internal reference gene were calculated by 2^−ΔΔCt (Yang et al. 2020). The primers targeting the representative DEGs were designed (Table S2), and synthesized by the Sangon (Shanghai, China).

Statistical analysis

The experimental data were analyzed using the SPSS software (version 17.0, SPSS Inc., Chicago, IL, USA). Differences between the means and changes in the samples were compared by one-way analysis of variance using the least-significant difference (LSD) method, with the level of significance set at p < 0.05. All tests were conducted in triplicate.

Results

Phenotypes of the V. parahaemolyticus isolates

V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 strains were used in this study, which were isolated from two species of shellfish (Solen strictus and Corbicula aurea), one species of crustacean (Oratosquilla oratoria), and one species of fish (Ctenopharyngodon idellus), respectively, in Shanghai, China in July and August of 2017 (Su and Chen 2020). The resistance of the strains to ten antibiotics were evaluated, and the results indicated that they were all resistant to ampicillin (AMP), kanamycin (KAN), rifampicin (RIF), and streptomycin (STR), while V. parahaemolyticus N2-5 was also resistant to tetracycline (TET) (Table 1). Meanwhile, these four isolates had different heavy metal-tolerance profiles (Table 1). For instance, V. parahaemolyticus B8-26 was tolerant to Hg²⁺/Cd²⁺/Cu²⁺/Zn²⁺, while V. parahaemolyticus B1-21 to Hg²⁺/Ni²⁺/Pb²⁺/Zn²⁺. In addition, V. parahaemolyticus ATCC17802 isolate of the human clinical origin was resistant to AMP (Jin et al. 2021), as well as Hg²⁺/Zn²⁺ (Melo-López et al. 2021).

Growth profiles of the V. parahaemolyticus isolates at different salinity concentrations, pH, temperatures, and oxygen concentrations

V. parahaemolyticus is a halophilic bacterium that inhabits in marine and estuarine environments worldwide (Li et al. 2019). We wondered growth features of the V. parahaemolyticus isolates of aquatic animal origins. Thus, we determined their growth curves in the TSB (pH 8.5) under the 0.5–5% NaCl conditions at 37 °C. As shown in Fig. S1 (A–E), all the V. parahaemolyticus isolates grew poorly at 0.5% NaCl. However, the increased NaCl concentration (1%) significantly promoted the bacterial growth (p < 0.05). Although all the isolates grew exuberantly at 2–5% NaCl, the highest biomass was observed when grew at 3% NaCl, showing the maximum OD₆₀₀ values (0.94–1.12) at stationary growth phase (SGP). Exceptionally, the growth of V. parahaemolyticus B1-21 was significantly decreased at 5% NaCl as compared to that at 3% NaCl (p < 0.05) (Fig. S1, B), suggesting that V. parahaemolyticus B1-21 was less tolerant to the higher NaCl concentration. These results indicated that the V. parahaemolyticus isolates grew optimally in the TSB at 3% NaCl.

The pH of human stomach normally ranges pH 1-3, but can rise above 6.0 after food consumption (Sun et al. 2014). The acidic stomach condition challenges the bacterial pass through to the gastrointestinal tract where it colonizes and causes the diarrhea disease. Thus, we determined growth curves of the V. parahaemolyticus isolates in the TSB (3% NaCl) under the pH 6.0–8.5 conditions at 37 °C. As shown in Fig. S2 (A–E), the growth of all the isolates was inhibited under the acidic conditions (pH 6.0–6.5), whereas the neutral condition (pH 7.0) promoted the bacterial growth, and all the isolates grew vigorously at pH 7.5–8.5, showing the maximum biomass at pH 8.5 (OD₆₀₀ = 0.94–1.20) at SGP. These results indicated that the V. parahaemolyticus isolates of aquatic animal origins grew optimally at pH 8.5, 3% NaCl in the TSB.

Growth curves of the V. parahaemolyticus isolates in the TSB (3% NaCl, pH 8.5) were also determined at 25 °C and 37 °C, respectively, wherein V. parahaemolyticus experiences during its life cycle in the environment and in the host. As shown in Fig. S3 (A–E), all the isolates could grow more exuberantly at 37 °C than at 25 °C. Therefore, the V. parahaemolyticus isolates were incubated in the TSB (pH 8.5, 3% NaCl) at 37 °C in the further analyses in this study, to avoid any condition influence other than the oxygen concentrations.

Growth of the V. parahaemolyticus isolates at different oxygen concentrations (21–5% O₂) was examined; and the results are presented in Fig. 1 A–E. Under the 5–10% O₂ conditions, the growth of all the isolates in the TSB (3% NaCl, pH 8.5) at 37 °C was repressed, and the maximum OD₆₀₀ values ranged 0.40–0.52 at SGP. Under the 12% O₂ condition, all the isolates grew better than at 10% O₂. Upon the increased O₂ concentrations (15–21%), the V. parahaemolyticus isolates grew faster correspondingly, but reaching the maximum biomass (OD₆₀₀ = 0.90–1.28) under the normal 21% O₂ condition. Additionally, we observed that the growth of V. parahaemolyticus ATCC17802 of the clinical origin was the slowest under the lower O₂ conditions among the test strains, moreover, it still showed the lower biomass under the 18% O₂ condition.

Fig. 1

The growth of the V. parahaemolyticus isolates in the TSB (pH 8.5, 3% NaCl) under different oxygen conditions at 37 °C. A–E: V. parahaemolyticus B8-26, V. parahaemolyticus B1-21, V. parahaemolyticus N2-5, V. parahaemolyticus L7-40, and V. parahaemolyticus ATCC17802, respectively

Effects of the hypoxic condition on biofilm formation of the V. parahaemolyticus isolates.

The ability of Vibrio to form biofilm attributes to their survival in the host and withstanding in different aquatic environments (Khan et al. 2020). Therefore, we asked whether and how the hypoxic condition (10% O₂) would affect biofilm formation of the V. parahaemolyticus isolates. The dynamic process of biofilm formation was monitored under the 10% O₂ condition when the isolates were incubated in the TSB (pH 8.5, 3% NaCl) at 37 °C for 60 h. As shown in Fig. 2A–E, the V. parahaemolyticus isolates formed biofilms at three different stages, including the development, maturation, and diffusion, similar to those under the normal 21% O₂ condition. However, different biofilm formation patterns were observed among the isolates.

Fig. 2

The biofilm formation of the V. parahaemolyticus isolates grown in the TSB (pH 8.5, 3% NaCl) under the 21% O₂ and 10% O₂ conditions at 37 °C. A–E: V. parahaemolyticus B8-26, V. parahaemolyticus B1-21, V. parahaemolyticus N2-5, V. parahaemolyticus L7-40, and V. parahaemolyticus ATCC17802, respectively. *p < 0.05; **p < 0.01; ***p < 0.001

Unexpectedly, under the 10% O₂ condition, the biofilm biomass generated by all the V. parahaemolyticus isolates was significantly increased at the development and maturation stages (0–36 h) as compared to those under the 21% O₂ condition (p < 0.05). The highest increase (2.70–fold) was observed in the V. parahaemolyticus B8-26 treatment group at 36 h (p < 0.001). Of note, the 10% O₂ condition significantly enhanced the biofilm formation of V. parahaemolyticus ATCC17802 for 48 h (p < 0.05) (Fig. 2, E). These results demonstrated that the hypoxic condition can enhance the biofilm formation of all the test V. parahaemolyticus isolates.

Cell morphological structure changes of the V. parahaemolyticus isolates under the hypoxic condition

Cell surface structure of the V. parahaemolyticus isolates under the hypoxic condition was further observed by the SEM assay. As shown in Fig. 3A–E, unexpectedly, in remarkable contrast to the control groups under the 21% O₂ condition, whose cell surface structure was broken with obvious crumpling, rupture and cellular contents leakage after 24-h incubation, the V. parahaemolyticus isolates showed rod cells with slightly shrinkage on the cell surface after being treated with 10% O₂ for 24 h. Particularly, V. parahaemolyticus ATCC17802 cells in the control group were severely shrunk and deformed, whereas those in the hypoxia treatment group were rod-shaped and less destroyed in cell structure (Fig. 3, E). These results indicated that the 10% O₂ condition repressed the growth of the V. parahaemolyticus isolates. The cellular metabolisms-requiring O₂ were consequently down-regulated, thereby leading to the slower growth under the hypoxic condition.

Fig. 3

The cell morphological structure of the V. parahaemolyticus isolates grown in the TSB (pH 8.5, 3% NaCl) under the 21% O₂ and 10% O₂ conditions at 37 °C for 24 h. A–E: V. parahaemolyticus B1-21, V. parahaemolyticus ATCC17802, V. parahaemolyticus B8-26, V. parahaemolyticus N2-5, and V. parahaemolyticus L7-40, respectively

Genome features of the V. parahaemolyticus isolates originating in edible aquatic animals

Based on the above results, draft genome sequences of V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 isolates of aquatic animal origins were determined, which yielded 37,042–51,777 clean sequencing reads. The assembled genome sizes were 4,913,675–5,353,490 bp with the GC contents of 45.29–45.44%. The yielded DNA scaffolds were 41–106. A total of 4,575–5,045 genes were predicted, of which the unknown function genes were 19.32–23.13% (Table 2, Fig. S4).

Table 2 Genome features of the V. parahaemolyticus isolates used in this study

V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 genomes also carried putative mobile genetic elements (MGEs), such as genome islands (GIs, n = 36), prophages (n = 4), integrons (INs, n = 26), and insertion sequences (ISs, n = 12), suggesting possible horizontal gene transfer (HGT) during the V. parahaemolyticus genome evolution. The draft genomes of V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 isolates were deposited in the GenBank database under the accession numbers JAODUR000000000, JAOPTY000000000, JAODVT000000000, and JAODVU000000000, respectively.

Additionally, based on the draft genome of V. parahaemolyticus ATCC17802 in the GenBank database (accession numbers NZ_CP014046.2) (Yang et al. 2015), we also identified several MGEs, including GIs (n = 2), and INs (n = 7) in this isolate.

GIs

GIs can shape bacterial genomes and affect bacterial fitness to the environment (Ghazali et al. 2021). In this study, 36 GIs were identified in the four V. parahaemolyticus genomes of aquatic animal origins, which contained 3–14 GIs with 5,572–38,679 bp, and carried 7–45 predicted genes with diverse biological functions (Fig. 4, Table S3).

Fig. 4

Gene organizations of the GIs identified in the V. parahaemolyticus genomes. Different colors referred to COG functional classification and genes with unknown function were in grey

The V. parahaemolyticus N2-5 genome contained the maximum number of GIs (n = 14, GIs 1–14), which were 5,572–38,427 bp, with 7–33 predicted genes. Conversely, the V. parahaemolyticus B8-26 genome carried the fewest GIs (n = 3, GIs 1–3), ranging 13,548–38,679 bp with 19–34 predicted genes. Notably, there were some GIs carrying virulence-related genes. For example, the GI 3 in V. parahaemolyticus B8-26, and the GI 11 in V. parahaemolyticus N2-5, encoding type II toxin-antitoxin (TA) system prophylactic-host death family antitoxins (Vp B8-26_3404, Vp N2-5_3709). The TA systems are involved in biofilm formation, persistence, stress endurance, and programmed cell death in favor of bacterial population (Hosseini et al. 2019).

Prophages

Phages, being the most abundant biological entities on earth, are viruses that infect bacteria, and constantly reshape bacterial communities (Wahl et al. 2019). In this study, we found four prophages in the V. parahaemolyticus B1-21, N2-5, and L7-40 genomes (Table S4), ranging 11,738–32,968 bp with 14–47 genes (Fig. 5). Like V. parahaemolyticus B8-26, no prophage was reported in V. parahaemolyticus ATCC17802 (Yang et al. 2015).

Fig. 5

The structure diagram of the prophages identified in the V. parahaemolyticus genomes

V. parahaemolyticus B1-21 genome contained two prophage gene clusters that showed high sequence similarity to Vibrio_phage_K139 (33,106 bp, NCBI accession number: NC_003313), and Vibrio_phage_fs2 (8,651 bp, NCBI accession number: NC_001956), respectively. The Vibrio_phage_K139 homologue was also found in the V. parahaemolyticus L7-40 genome, but with a truncated version (18,896 bp) carrying 23 genes, suggesting the rearrangement of Vibrio_phage_K139 during the V. parahaemolyticus genome evolution. The V. parahaemolyticus N2-5 genome contained a truncated version (29,468 bp) of Enterobacteria_phage_Mu (36,717 bp, NCBI accession number: NC_000929), carrying 47 genes encoding 32 unknown proteins, which suggested possible transmission of the Enterobacteria_phage_Mu between the Vibrio and Enterobacteria genera pathogens.

INs

INs are also critically involved in bacterial evolution and antimicrobial resistance with the carried gene cassettes (Ali et al. 2022). In this study, all the V. parahaemolyticus genomes contained INs (n = 2–11) ranging 605 bp–65,657 bp (Fig. S5-Fig. S6, Table S5).

The V. parahaemolyticus B1-21 genome contained the maximum number of INs (n = 11, INs 1–11), but all of which were incomplete INs with gene cassettes. In contrast, the V. parahaemolyticus L7-40 genome carried only two INs (INs 1–2), including one complete IN (IN 1) and one gene cassette (IN 2). Additionally, seven INs were identified in the V. parahaemolyticus ATCC17802 genome, including one complete (IN 1) with an integrase IntI 4 (Vp B8-26_1368), suggesting that the IN 1 was a super integron (Yang et al. 2015).

ISs

ISs are mobile repeat sequences that can copy themselves to new locations in bacterial genomes (Tempel et al. 2022). In this study, ISs were identified in the V. parahaemolyticus B8-26 (n = 4), B1-21 (n = 1), N2-5 (n = 5), and L7-40 (n = 2) genomes (Table S6), but absent in V. parahaemolyticus ATCC17802 (Yang et al. 2015). Notably, one IS001 (1,307 bp) belonging to the ISAs1 family was identified in V. parahaemolyticus L7-40, encoding a c-di-GMP metabolism protein (Vp L7-40_4574), which is critical to bacterial adaptation to changing environments, including biofilm formation, host colonization and virulence (Kumar and Chatterji 2008).

Putative virulence and resistance-associated genes

Many virulence-related genes (n = 42–45) were identified in the V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 genomes. For instance, more than thirty genes for type III secretion system 1 (T3SS1) were found, encoding the VecA, TyeA, SycN, YscO, VcrDGHRV, VopBDNQRS, and VscBCDFGHIJKLNQRSTUXY. The ExsA gene was also found in these genomes (Vp B8-26_3495, Vp B1-21_3659, Vp N2-5_3801, and Vp L7-40_3386), which is the master transcription factor that positively regulates T3SS1 expression (Gu et al. 2020). In addition, V. parahaemolyticus ATCC17802 has both T3SS1 and T3SS2, and TRHs (Yang et al. 2015).

Several antibiotic resistance-related genes (n = 3) and heavy metal tolerance-ralated genes (n = 8) were identified in the V. parahaemolyticus B8-26, B1-21, N2-5, L7-40, and ATCC17802 genomes (Table 3), which provided the genome-wide evidence for resistance phenotypes of the V. parahaemolyticus isolates.

Table 3 The antibiotic and heavy metal resistance-related genes identified in the V. parahaemolyticus genomes

Comparative transcriptome analysis of the V. parahaemolyticus isolates in response to the hypoxic condition

Based on the V. parahaemolyticus B8-26, B1-21, N2-5, L7-40 genomes determined in this study, coupled with the V. parahaemolyticus ATCC17802 genome, we further examined genome-wide gene expression changes triggered by the hypoxic condition. The isolates were incubated in the TSB under the 10% O₂ condition at 37 °C for 24 h, at which the cell morphological structure changes of the V. parahaemolyticus isolates was the most obvious, therefore, their transcriptomes were determined. The complete lists of DEGs in the V. parahaemolyticus isolates are available in the NCBI SRA database (https://submit.ncbi.nlm.nih.gov/subs/bioproject/) under the accession number PRJNA906699 and PRJNA767551.

The major altered metabolic pathways in V. parahaemolyticus B8-26 triggered by the hypoxic condition

Approximately 65.89% (2,722/4,131) of V. parahaemolyticus B8-26 genes were expressed differently under the hypoxic condition as compared to the control group. Among these, 413 differently expressed genes (DEGs) showed higher transcription levels (fold change, FC ≥ 2.0), whereas 2,309 DEGs were down-regulated (FC ≤ 0.5). The comparative transcriptomic analysis revealed fourteen significantly changed metabolic pathways, including the aminoacyl-tRNA biosynthesis; protein export; peptidoglycan (PG) biosynthesis; citrate cycle; sulfur metabolism; ubiquinone and other terpenoid-quinone biosynthesis;, RNA degradation; glyoxylate and dicarboxylate metabolism; purine metabolism; thiamine metabolism; ATP-binding cassette (ABC) transporters; pyruvate metabolism; tyrosine metabolism; and carbon fixation pathways in prokaryotes (Fig. 6, Table S7).

Fig. 6

The Volcano plot of differential gene expression (A), and the major changed metabolic pathways (B) in V. parahaemolyticus strains under the hypoxic condition (10% O₂). A: The X and Y axes represented changes of the up/down-regulated genes, and the corresponding significant differences, respectively. B The X and Y axes represented the major changed metabolic pathways in V. parahaemolyticus strains under the hypoxic condition, and the ratios of the number of DEGs in the metabolic pathways to the total number of DEGs in all metabolic pathways. A1 and B1, A2 and B2, A3 and B3, A4 and B4: V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 strains, respectively

Remarkably, all DEGs in the aminoacyl-tRNA biosynthesis (n = 25), PG biosynthesis (n = 17), citrate cycle (n = 8), and ubiquinone and other terpenoid-quinone biosynthesis (n = 17) were significantly repressed at the transcription levels in V. parahaemolyticus B8-26 under the hypoxic condition (0.015–fold to 0.497–fold) (p < 0.05). Aminoacyl-tRNA synthetases are essential for protein synthesis (Rubio Gomez and Ibba 2020). The expression of many such enzymes and proteins was repressed in V. parahaemolyticus B8-26. For instance, the DEG encoding an alanyl-tRNA editing protein (Vp_B8-26_03025) was strongly inhibited (0.080–fold), suggesting that the stringent response shut down the translation to avoid toxic generation of mistranslated/misfolded proteins (Aggarwal et al. 2021). As the major component of bacterial cell wall, the biosynthesis of PG is closely related to cell wall growth and turnover (Wamp et al. 2020). In this study, seventeen DEGs involved in the PG biosynthesis were inhibited under the hypoxic condition. For example, the DEG encoding a D-alanyl-D-alanine carboxypeptidase (Vp_B8-26_24445) was significantly down-regulated (0.448–fold), which is involved in the reconstruction of newly synthesized PG in Francisella (Spidlova et al. 2018). The citrate cycle (Krebs cycle) is a major aerobic pathway for the final stages of carbohydrate and fatty acid oxidation, which provides NADH to be used for oxidative phosphorylation and other metabolic reactions. In this study, eight DEGs in the citrate cycle were down-regulated as well. For instance, the DEG encoding a malate dehydrogenases (MDH) (Vp_B8-26_01555) was strongly down-regulated (0.022–fold). The MDH catalyzes the NAD/NADH-dependent interconversion of malate and oxaloacetate. The malate/aspartate of this reaction is involved in tricarboxylic acid cycle (Minárik et al. 2002). Additionally, among the seventeen down-regulated DEGs in the ubiquinone and other terpenoid-quinone biosynthesis, the DEG encoding a NAD(P)H quinone oxidoreductase (NQO) (Vp_B8-26_12250) was significantly repressed (0.172–fold), which is a two-electron reductase responsible for detoxification of quinones and bioactivation of certain quinones (Zhang et al. 2018). These results signified that the hypoxic condition likely inhibited protein translation, bacterial cell wall synthesis, and energy supply in V. parahaemolyticus B8-26.

Sulfur is also a requirement for both the host and colonizing bacteria (Kies et al. 2022). In this study, fourteen DEGs were significantly down-regulated in the sulfur metabolism in V. parahaemolyticus B8-26 under the hypoxic condition (0.010–fold to 6.247–fold) (p < 0.05). For example, the NrfD-like subunits are found in many diverse membrane complexes, which may participate in the metabolism of oxygen, nitrogen, sulfur, arsenate or hydrogen (Calisto and Pereira 2021). Nevertheless, in this study, the expression of NrfD (Vp_B8-26_10960) was strongly down-regulated (0.007–fold), suggesting the inhibited sulfur metabolism, thereby energy conversion in V. parahaemolyticus B8-26 under the hypoxic condition.

Of note, approximately 85 DEGs in the ABC transport were also significantly inhibited in V. parahaemolyticus B8-26 under the hypoxic condition (0.003–fold to 0.498–fold) (p < 0.05). The ATP-dependent ABC transporter proteins facilitate the import and/or export of various substrates, including lipids, sugars, amino acids and peptides, ions, and drugs (Thurm et al. 2021). For example, the DEG encoding the ATP-binding cassette domain-containing protein (Vp_B8-26_06880) was greatly down-regulated (0.003–fold), which contains two domains capable of ATP hydrolysis in order to fuel protein function (Alqahtani et al. 2019). The down-regulation of the ABC transport suggested the inactive cellular transport function under the hypoxia, which likely led to harmful substance accumulate and retarded bacterial growth.

Conversely, the expression of some DEGs was significantly increased in V. parahaemolyticus B8-26 under the hypoxic condition (p < 0.05). For instance, four DEGs involved in the protein export were significantly up-regulated (2.201–fold to 7.248–fold) (p < 0.05), e.g., the protein translocase subunit SecD (Vp_B8-26_23120, 7.248-fold), and protein translocase subunit SecF (Vp_B8-26_23115, 5.060-fold). The SecDF belong to the SecYEG-SecDF-YajC-YidC holo-translocon (HTL) protein secretase/insertase, a supercomplex necessary for protein secretion, insertion membrane and complex assembly (Komar et al. 2016).

Taken, the comparative transcriptome data revealed fourteen significantly changed metabolic pathways in V. parahaemolyticus B8-26 under the hypoxic condition. Notably, the down-regulated protein translation and transport, ABC transporters, bacterial cell wall synthesis and energy conversion likely contributed to the retarded cell growth of V. parahaemolyticus B8-26 under the hypoxic condition.

The major changed metabolic pathways in V. parahaemolyticus B1-21 triggered by the hypoxic condition

Approximately 46.77% (2,058/4,400) of V. parahaemolyticus B1-21 genes were expressed differently under the hypoxic condition as compared to the control group. Of these, 1,533 DEGs were up-regulated (FC ≥ 2.0), whereas 525 DEGs were down-regulated (FC ≤ 0.5). Eleven significantly changed metabolic pathways were identified, including the ribosome; phosphotransferase system (PTS); two-component system (TCS); bacterial chemotaxis; ABC transporters; alanine, aspartate and glutamate metabolism; flagellar assembly; plant-pathogen interaction; sulfur relay system; sulfur metabolism; and pyruvate metabolism (Fig. 6, Table S8).

Bacterial chemotaxis is a critical ability to search for the optimal environment to ensure the survival of bacterial species (Jeong 2021). Three DEGs in the bacterial chemotaxis were significantly repressed in V. parahaemolyticus B1-21 under the hypoxic condition (0.188–fold to 0.465–fold) (p < 0.05), including a chemotactic protein Chew (Vp_B1-21_11985), a flagellar motor switch protein (Vp_B1-21_12100), and a chemotactic response regulator protein-glutamate methylase (Vp_B1-21_12000). For example, the DEG encoding the CheW (Vp_B1-21_11985), which is involved in the transmission of sensory signals from the chemoreceptors to the flagellar motors (Rosario et al. 1994), was significantly repressed (0.465–fold). Meanwhile, the DEG encoding the flagellar motor switch protein (Vp_B1-21_12100) was also highly repressed (0.188–fold). These results indicated the down-regulated bacterial chemotaxis under the hypoxic condition, thereby inhibiting bacterial migration to the environment conducive to survival.

Approximately 28 DEGs in the TCS were significantly inhibited (0.053–fold to 0.494–fold), whereas 56 DEGs were significantly enhanced (2.017–fold to 0.494–fold) at the transcriptional levels in V. parahaemolyticus B1-21 under the hypoxic condition (p < 0.05). Bacteria sense and respond to environmental changes via TCS, e.g., cell surface modifications, and increased biofilm formation (Tierney and Rather 2019). For example, the DEG encoding an acetyl coenzyme A acetyltransferase (ACAT) (Vp_B1-21_23095, 0.121–fold) was significantly inhibited in V. parahaemolyticus B1-21 (p < 0.05), which catalyzes the reversible formation of acetylacyl coenzyme A from two acetyl coenzyme A molecules during the ketogenesis and ketolysis (Goudarzi 2019).

Remarkably, all DEGs (n = 42) in the ribosome were significantly down-regulated (0.214–fold to 0.495–fold) in V. parahaemolyticus B1-21 under the hypoxic condition (p < 0.05), e.g., 50S ribosomal proteins (Vp_B1-21_01230, Vp_B1-21_06585, Vp_B1-21_13890), and 30S ribosomal proteins (Vp_B1-21_01360, Vp_B1-21_11045, Vp_B1-21_15475, Vp_B1-21_15315). Ribosomes are essential for protein production, and thus for bacterial growth, and proliferation (Turi et al. 2019). These results indicated that V. parahaemolyticus B1-21 reduced protein translation in response to the hypoxic condition.

Conversely, all DEGs (n = 21) in the PTS were significantly up-regulated (2.616–fold to 28.083–fold) in V. parahaemolyticus B1-21 under the hypoxic condition (p < 0.05). For instance, the gene encoding a fused PTS fructose transporter subunit IIA/HPr protein (Vp_B1-21_21185) was highly up-regulated (28.083–fold), suggesting the up-regulated transport and utilization of carbohydrates in V. parahaemolyticus B1-21 under the hypoxic condition.

Similar to V. parahaemolyticus B8-26, the ABC transporters, sulfur metabolism, and pyruvate metabolism were also significantly changed in V. parahaemolyticus B1-21 under the hypoxic condition.

Taken together, the comparative transcriptome data revealed eleven significantly changed metabolic pathways in V. parahaemolyticus B1-21 under the hypoxic condition. Specially, V. parahaemolyticus B1-21 reduced protein translation, but enhanced transport and utilization of carbohydrates in response to the hypoxic condition.

The major changed metabolic pathways in V. parahaemolyticus N2-5 triggered by the hypoxic condition

Approximately 47.55% (783/4,400) of V. parahaemolyticus N2-5 genes were expressed differently under the hypoxic condition as compared to the control group. Of these, 1,102 DEGs were up-regulated (FC ≥ 2.0), whereas 965 DEGs were down-regulated (FC ≤ 0.5). Eleven significantly changed metabolic pathways were identified, including ABC transporters; the ribosome; butanoate metabolism; tryptophan metabolism; alanine, aspartate and glutamate metabolism; propanoate metabolism; pyruvate metabolism; oxidative phosphorylation; valine, leucine and isoleucine degradation; citrate cycle; and carbon fixation pathways in prokaryotes (Fig. 6, Table S9).

The same case as V. parahaemolyticus B1-21, all DEGs (n = 45) in the ribosomes were significantly down-regulated (0.137–fold to 0.468–fold) at the transcription levels in V. parahaemolyticus N2-5 under the hypoxic condition (p < 0.05). Like V. parahaemolyticus B8-26, all DEGs (n = 5) in the citrate cycle were also significantly repressed (0.040–fold to 0.188–fold) in V. parahaemolyticus N2-5.

Of note, almost all DEGs (n = 24) in the oxidative phosphorylation were significantly repressed (0.062–fold to 0.497–fold) in V. parahaemolyticus N2-5 under the hypoxic condition (p < 0.05), e.g., the F0F1 ATP synthase subunits A, B, α, β, and γ (Vp_N2-5_16940, Vp_N2-5_16930, Vp_N2-5_16920, Vp_N2-5_16910, and Vp_N2-5_16915), a cytochrome o ubiquinol oxidase subunit I (Vp_N2-5_20305), a cytochrome d ubiquinol oxidase subunit II (Vp_N2-5_05315), a cytochrome o ubiquinol oxidase subunit III (Vp_N2-5_20310), a cytochrome ubiquinone oxidase (Vp_N2-5_05315), a cytochrome c1, a cytochrome bc complex cytochrome b subunit (Vp_N2-5_02185), and a protoheme IX farnesyltransferase (Vp_N2-5_20320). The oxidative phosphorylation provides most of the ATP that is required for setting and maintaining cellular metabolic homeostasis (Wilson 2017). The significant repressed oxidative phosphorylationof indicated insufficient cellular energy production in V. parahaemolyticus N2-5, consequently affecting energy expenditure and cell growth under the hypoxic condition.

In the butanoate metabolism, thirteen DEGs were significantly down-regulated (0.028–fold to 0.475–fold) in V. parahaemolyticus N2-5 (p < 0.05), whereas five DEGs were significantly up-regulated (3.077–fold to 14.428–fold) (p < 0.05). For example, the DEG encoding a class I poly(R)-hydroxyalkanoic acid (PHA) synthase (Vp_N2-5_23085) was highly down-regulated (0.028–fold). PHAs are aliphatic polyesters produced by many bacteria and archaea in response to various environmental conditions (McCool and Cannon 2001).

In the tryptophan metabolism, three DEGs were significantly down-regulated (0.033–fold to 0.347–fold) in V. parahaemolyticus N2-5 (p < 0.05), whereas two DEGs were significantly up-regulated (8.781–to 10.778–fold) (p < 0.05). Of these, the expression of a catalase (Vp_N2-5_23615) was highly up-regulated (8.781–fold). The catalase decomposes hydrogen peroxide, protecting cells from potentially harmful reactive oxygen species (Kim et al. 2018).

Similar to V. parahaemolyticus B1-21, the ABC transport, and alanine, aspartate and glutamate metabolism were significantly changed in V. parahaemolyticus N2-5 as well.

Taken, the comparative transcriptome analysis revealed eleven significantly changed metabolic pathways in V. parahaemolyticus N2-5 under the hypoxic condition, some of which were also altered in V. parahaemolyticus B1-21 and V. parahaemolyticus B8-26. Specially, the protein translation, citrate cycle, propionate metabolism, oxidative phosphorylation, carbon fixation pathways in prokaryotes, and valine, leucine and isoleucine degradation were all repressed in V. parahaemolyticus N2-5 for growth under the hypoxic condition.

The major changed metabolic pathways in V. parahaemolyticus L7-40 triggered by the hypoxic condition

Approximately 37.43% (783/4,400) of V. parahaemolyticus L7-40 genes were expressed differently under the hypoxic condition as compared to the control group. Of these, 1,032 DEGs were up-regulated (FC ≥ 2.0), whereas 615 DEGs were repressed (FC ≤ 0.5). Fifteen significantly changed metabolic pathways were identified, including the arginine biosynthesis; glyoxylate and dicarboxylate metabolism; glycolysis/gluconeogenesis; methane metabolism; alanine, aspartate and glutamate metabolism; pyruvate metabolism; oxidative phosphorylation; propanoate metabolism; citrate cycle; beta-lactam resistance; TCS; taurine and hypotaurine metabolism; PTS; nitrogen metabolism; and carbon fixation pathways in prokaryotes (Fig. 6, Table S10).

Like V. parahaemolyticus B8-26, and N2-5 strains, all DEGs (n = 5) in the citrate cycle were significantly down-regulated (0.106–fold to 0.471–fold) in V. parahaemolyticus L7-40 under the hypoxic condition (p < 0.05). For example, the DEG encoding a fumarate reductase (Vp_L7-40_15840) was significantly inhibited (0.471–fold) (p < 0.05), which is a flavin protease containing the cofactor flavin adenine dinucleotide and reduces fumarate to succinate (Alqahtani et al. 2019). Fumarate reductases are essential for maintaining redox homeostasis in cells because they reoxidize intracellular flavin adenine dinucleotides (FADH₂) under the hypoxic conditions (Kim et al. 2018). The DEG encoding a succinate dehydrogenase (SDH) (Vp_L7-40_15845) was also inhibited (0.111–fold) (p < 0.05), which is a critical enzyme involved in the tricarboxylic acid cycle and oxidative phosphorylation for energy production (Moosavi et al. 2019). These results suggested that the hypoxia likely disrupted the cellular redox balance and the energy supply in V. parahaemolyticus L7-40.

Like V. parahaemolyticus N2-5, all DEGs (n = 14) in the propanoate metabolism were significantly down-regulated (0.121–fold to 0.437–fold) in V. parahaemolyticus L7-40 (p < 0.05). The same case as V. parahaemolyticus N2-5, the 16 of 21 DEGs in the oxidative phosphorylation were significantly down-regulated (0.116–fold to 0.454–fold) as well (p < 0.05). Interestingly, in marked contrast to the most inhibited DEGs in the oxidative phosphorylation, the expression of one DEG encoding a manganese-dependent inorganic pyrophosphatase (PPase, Vp_L7-40_09185), which plays an essential role in the conservation of energy and supplies energy for numerous biosynthetic pathways (Hu et al. 2020), was strongly up-regulated (33.120–fold) in V. parahaemolyticus L7-40 (p < 0.05).

In the glyoxylate and dicarboxylate metabolism, the 20 of 24 DEGs were significantly down-regulated (0.019–fold to 0.491–fold) in V. parahaemolyticus L7-40 under the hypoxic condition (p < 0.05). Of these, the DEG encoding an acetate-CoA ligase (Vp_L7-40_16020) was highly down-regulated (0.019–fold). Acetyl-CoA is a vitally important and versatile metabolite used for many cellular processes, and can also deal with stress such as low nutritional availability and hypoxia (Miller et al. 2021).

Similar to V. parahaemolyticus B1-21, all DEGs (n = 13) in the PTS were significantly up-regulated (2.338–fold to 64.585–fold) in V. parahaemolyticus L7-40 under the hypoxic condition (p < 0.05). In the glycolysis/gluconeogenesis, the 17 of 21 DEGs were significantly up-regulated (2.199–fold to 10.689–fold) in V. parahaemolyticus L7-40 (p < 0.05). Of these, the DEG encoding a 6-phosphofructokinase (PFK, Vp_L7-40_15910) was highly up-regulated (10.689–fold). The PFK is one of the most prominent rate-limiting enzymes in the glycolytic pathway. The increase in glycolytic flux is beneficial for maintaining bioenergetic homeostasis during the hypoxia (Kierans and Taylor 2021).

In addition, similar to V. parahaemolyticus B8-26, and V. parahaemolyticus N2-5, the TCS, pyruvate metabolism, as well as alanine, aspartate and glutamate metabolism were also significantly changed in V. parahaemolyticus L7-40 under the hypoxic condition.

Taken, the comparative transcriptome analysis revealed fifteen significantly changed metabolic pathways in V. parahaemolyticus L7-40, some of which were also identified in V. parahaemolyticus B8-26, N2-5, and B1-21 isolates under the hypoxic condition. Specifically, V. parahaemolyticus L7-40 inhibited the glyoxylate and dicarboxylate metabolism, but up-regulated the glycolysis/gluconeogenesis, and changed the nitrogen metabolism, and the arginine biosynthesis in response to the hypoxic condition.

The major changed metabolic pathways in V. parahaemolyticus ATCC17802 triggered by the hypoxic condition

Approximately 57.12% (783/4,674) of V. parahaemolyticus ATCC17802 genes were expressed differently under the hypoxic condition as compared to the control group. Among these, 361 DEGs were up-regulated (FC ≥ 2.0), whereas 2,309 DEGs were down-regulated (FC ≤ 0.5). Ten significantly changed metabolic pathways were identified, including the ribosome; ABC transporters; aminoacyl-tRNA biosynthesis; oxidative phosphorylation; pyrimidine metabolism; fatty acid biosynthesis; arginine and proline metabolism; lipopolysaccharide (LPS) biosynthesis; purine metabolism, and carbon fixation pathways in prokaryotes (Fig. 7, Table S11).

Fig. 7

The Volcano plot of differential gene expression (A), and the major changed metabolic pathways (B) in V. parahaemolyticus ATCC17802 under the hypoxic condition

Similar to V. parahaemolyticus B1-21, and N2-5 isolates, all DEGs (n = 49) in the ribosome were significantly down-regulated (0.033–fold to 0.243–fold) at the transcriptional level in V. parahaemolyticus ATCC17802 under the hypoxic condition (p < 0.05). The same case as V. parahaemolyticus N2-5, and L7-40 isolates, all DEGs (n = 35) in the oxidative phosphorylation were significantly down-regulated (0.029–fold to 0.440–fold) in V. parahaemolyticus ATCC17802, except the polyphosphate kinase (Vp_17802_02840, 2.187–fold).

In the fatty acid biosynthesis, all DEGs (n = 17) were significantly down-regulated (0.047–fold to 0.427–fold) in V. parahaemolyticus ATCC17802 (p < 0.05). Fatty acid biosynthesis plays a central role in building cell membrane, reserving cell energy, and production of precursors to second messenger molecules (Günenc et al. 2022). For instance, the DEG encoding a trans-2-enoyl-CoA reductase family protein (Vp_17802_21460) was strongly down-regulated (0.047–fold) in V. parahaemolyticus ATCC17802, which is involved in chain lengthening of fatty acids (Uchida et al. 2021).

In the LPS biosynthesis, most DEGs (16/18) were significantly down-regulated (0.068–fold to 0.476–fold) in V. parahaemolyticus ATCC17802 (p < 0.05). The LPS, the major component of the outer membrane of Gram-negative bacteria, is essential for bacterial viability and host–pathogen interactions (Kutschera et al. 2021). For instance, the DEG encoding the tetraacyldisaccharide 4′-kinase (Vp_17802_04990) was strongly down-regulated (0.068–fold) in V. parahaemolyticus ATCC17802, which is the prime enzyme for lipid biosynthesis vital for bacterial survival (Damale et al. 2022).

In the pyrimidine metabolism, most DEGs (26/30) were significantly down-regulated (0.026–fold to 0.450–fold) in V. parahaemolyticus ATCC17802 (p < 0.05). Pyrimidines are structural elements of a wide range of essential compounds in the synthesis of DNA, RNA, lipids, and carbohydrates (Garavito et al. 2015). For instance, the DEG encoding a thymidine kinase (TK, Vp_17802_05700) was strongly down-regulated (0.051–fold) in V. parahaemolyticus ATCC17802. The TK is involved in the pyrimidine nucleotide recovery pathway, a cell-proliferation marker (Fanelli et al. 2021). Interestingly, among the four up-regulated DEGs in the pyrimidine metabolism, the DEG encoding a 5' deoxynucleotidase (Vp_17802_04570, 4.317–fold) was highly up-regulated in V. parahaemolyticus ATCC17802 (p < 0.05). This enzyme is involved in the metabolism of cladribine (2CdA), by mediating phosphorolysis (deactivation) of 2CdA to prevent its accumulation in the cell, which leads to cell death (Carlini et al. 2022).

The same cases as V. parahaemolyticus B1-21, and N2-5 isolates, all the DEGs (n = 49) in the ribosome were highly down-regulated (0.033–fold to 0.243–fold) in V. parahaemolyticus ATCC17802 under the hypoxic condition (p < 0.05). Similar to V. parahaemolyticus B8-26, B1-21, and N2-5 isolates, almost all DEGs (105/108) in the ABC transporters were also significantly repressed (0.001–fold to 0.479–fold) in V. parahaemolyticus ATCC17802 (p < 0.05). Like V. parahaemolyticus B8-26, N2-5, and L7-40 isolates, most DEGs (19/21) in the carbon fixation pathways in prokaryotes were significantly inhibited (0.049–fold to 0.359–fold) in V. parahaemolyticus ATCC17802 (p < 0.05). In addition, like V. parahaemolyticus B8-26, all DEGs (n = 25) in the aminoacyl-tRNA biosynthesis were significantly repressed (0.100–fold to 0.477–fold) in V. parahaemolyticus ATCC17802 as well (p < 0.05).

Taken, the comparative transcriptome analysis revealed ten significantly changed metabolic pathways in V. parahaemolyticus ATCC17802 under the hypoxic condition. Although some of the altered pathways were also observed in V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 isolates, however, wherein more repressed DEGs were elicited in V. parahaemolyticus ATCC1780 by the condition, consistent with the phenotypes of this isolate in the above results. In particular, compared to the other four V. parahaemolyticus isolates, V. parahaemolyticus ATCC17802 also down-regulated some metabolic pathways, including the pyrimidine metabolism, fatty acid biosynthesis, lipopolysaccharide biosynthesis, purine metabolism, as well as arginine and proline metabolism to survival under the hypoxic condition.

Additionally, to confirm the transcriptome data, thirty representative DEGs were analyzed by the RT-qPCR assay, and the resulting data were generally correlated with those by the transcriptome analysis (Table S12).

Possible molecular mechanisms of V. parahaemolyticus in response to the hypoxic condition

The transcriptome-wide analyses revealed a number of DEGs involved in multiple pathways of biosynthesis, degradation, salvage, interconversion, and transport of the compounds in V. parahaemolyticus B8-26, B1-21, N2-5, L7-40 and ATCC17802 isolates of aquatic animal and human clinical origins, suggesting a complex molecular regulation network in the bacterium in response to the hypoxic condition (Fig. S7, Table S13).

Given that the V. parahaemolyticus isolates, being from different origins, harbored different genome features, it was reasonable that the hypoxic condition triggered different transcriptomic profiles among the isolates. The transcriptome data, coupled with the phenotype results in this study, demonstrated that compared to the other isolates, V. parahaemolyticus ATCC17802 down-regulated and/or shut down all the ten changed metabolic pathways to reduce cell viability and maintain cell structure integrated under the hypoxic condition. These results provided a novel mechanism for the persistence of the toxic strain in the environment and in the host.

On the other hand, same identical metabolic pathways were elicited in the V. parahaemolyticus isolates by the hypoxic condition, e.g., the inhibited ribosome in V. parahaemolyticus B1-21, N2-5, and ATCC17802 isolates; the repressed aminoacyl-tRNA biosynthesis in V. parahaemolyticus B8-26, and ATCC17802 isolates; and the repressed oxidative phosphorylation in V. parahaemolyticus N2-5, L7-40, and ATCC17802 isolates.

Overall, V. parahaemolyticus developed multiple molecular strategies to efficiently mitigate the cell damage and/or cytotoxicity caused by the hypoxia: (1) down-regulated the oxidative phosphorylation, carbon fixation pathways, citrate cycle, pyruvate metabolism, and propionate metabolism to balance the redox state of the cell and the energy conversion; (2) down-regulated the fatty acid biosynthesis, and lipopolysaccharide biosynthesis to inhibit the cell wall synthesis, thereby the cell proliferation; (3) decreased amino acid transport, aminoacyl-tRNA biosynthesis, and protein production to retard the cell growth, thus prolonging the growth cycle and maintaining the cell structure integrated; and (4) conversely, enhanced the expression of condition-related proteins (such as GapA), structurally stabilizing factors (such as arginine), and efflux RND transporters to reduce the cell damage for growth under the unfavorable hypoxic condition.

Serotypes of the V. parahaemolyticus isolates

The BLAST analysis of the antigen gene loci reveled that the V. parahaemolyticus B8-26 and B1-21 genomes contained the O antigen loci wvaG/wvcA, and specific loci VP6 and VP202 for K34 and K23 polymorphic sites, respectively, indicating that the serotypes of V. parahaemolyticus B8-26, and B1-21 isolates were O1:K34, and O5:K23, respectively. Similarly, V. parahaemolyticus N2-5, and L7-40 genomes carried O antigen loci orf16/wzy and wvaH, indicating that the serotypes of V. parahaemolyticus N2-5, and L7-40 isolates were O4/O11:K4, and O9:KUT, respectively. Additionally, V. parahaemolyticus ATCC17802 carried O antigen loci gmd, and specific loci VP199 for K1 polymorphic sites, and its serotype was O1:K1 (Bian et al. 2020).

Phylogenetic relatedness of the V. parahaemolyticus isolates

To address the phylogenetic relatedness of the V. parahaemolyticus isolates, we constructed a phylogenetic tree based on a total of 1,921 homologous amino acid sequences identified from the seventy-eight V. parahaemolyticus genomes, seventy-four of which were derived from the GenBank database (Table S1). This analysis revealed ten different phylogenetic groups, designated as Groups 1–10 (Fig. 8).

Fig. 8

The phylogenetic tree showing the relationship of the seventy-eight V. parahaemolyticus genomes. Draft genome sequences of the V. parahaemolyticus isolates determined in this study were marked with red dots in the tree

V. parahaemolyticus B8-26 (serotype: O1:K34; GenBank accession number: JAODUR000000000) was classified into Group 8, showing the closest evolutionary distance with V. parahaemolyticus XMM117 (serotype: O3:K6; GenBank accession number: NZ_CP064037.1), which was isolated from the environment in 2019 in China.

V. parahaemolyticus L7-40 (serotype: O9:KUT; GenBank accession number JAODVU000000000) fell into Group 4, together with V. parahaemolyticus FDAARGOS-667 (serotype: O5:K33/K55; GenBank accession number: NZ_CP044062.1), which was isolated from Homo sapiens in USA (collection time unknown).

V. parahaemolyticus B1-21 (serotype: O5:K23; GenBank accession number: JAOPTY000000000), and V. parahaemolyticus N2-5 (serotype: O4/O11:K4; GenBank accession number: JAODVT000000000) were classified into single evolutionary branches Groups 3, and 9, respectively, indicating that these two isolates have unique genome features.

In addition, the V. parahaemolyticus strains that can cause gastroenteritis in humans were classified into Groups 1, 2, 5, 6, 7, and 8, e.g., V. parahaemolyticus ATCC17802 (serotype: O1:K1; GenBank accession number: NZ_CP014046.1) in Group 2; V. parahaemolyticus 64 (serotype: O1:KUT; GenBank accession number: NZ_CP074415.1) isolated from Penaeus in Group 1; and V. parahaemolyticus XMM117 (serotype: O3:K6; GenBank accession number: NZ_CP064037.1) isolated from the environment in Group 8.

Taken, these results demonstrated the diversity of V. parahaemolyticus genomes of the environmental and clinical origins.

Discussion

Infectious diseases caused by pathogenic bacteria continue to be a global concern for public health, which results in millions of deaths worldwide each year (Bueno et al. 2020). The sea foodborne illness caused by V. parahaemolyticus has shown a significant upward trend in recent years (Zhai et al. 2021). Climate change is driving ocean oxygen declines (Jaccard and Galbraith 2012). On the other hand, the low oxygen level in the human gastrointestinal tract also challenges the survival and infection of V. parahaemolyticus. Thus, in the present study, we for the first time investigated how V. parahaemolyticus isolates of aquatic animal and human clinical origins response to the hypoxic condition.

V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 strains were isolated from aquatic animals S. strictus, C. aurea, O. oratoria, and C. idellus, respectively (Su and Chen 2020). Multiple antibiotic and heavy metal resistance profiles were derived from the four V. parahaemolyticus isolates, suggested possible antibiotic and heavy metal exposure or pollution sources in the aquaculture environments (Su and Chen 2020; Xu et al. 2022a). Moreover, we observed that these isolates grew optimally in the TSB medium with 3% NaCl and pH 8.5 at 37 °C, consistent with the previous report (Xu et al. 2022a).

V. cholerae is a member of the Vibrio genus and can cause pandemic cholerae in humans (Ramamurthy et al. 2020). V. cholerae was able to generate energy and maintain its physiological functions in the absence of oxygen using alternative electron acceptors (AEAs) (Bueno et al. 2020). In this study, we determined growth curves of the V. parahaemolyticus isolates in the TSB equalibrated with different oxygen concentrations (21–5% O₂) at 37 °C. Unlike V. cholerae, our data showed that the lower oxygen concentrations (5–10%) greatly inhibited the growth of the tested isolates. However, no significant difference in the growth was observed between 5% O₂ and 10% O₂, suggesting that 10% O₂ was the cut-off point for the growth of V. parahaemolyticus at the low oxygen concentrations. Additionally, we observed that oxygen concentrations in air within the culture tubes fluctuated down by 0–1% within 24 h, and by 0–2% from 24 to 48 h, indicating the stable hypoxic condition maintained in our culture system.

Biofilm is critical for V. parahaemolyticus persistence in aquatic environments and pathogenicity in the host (Yildiz and Visick 2009). In this study, the results showed that the biofilm biomass generated by all the tested V. parahaemolyticus isolates under the 10% O₂ condition was significantly increased at the development and maturation stages (0–36 h) (p < 0.05), as compared to those under the normal 21% O₂ condition. It has been reported that biofilm formation of P. aeruginosa, S. aureus, and E. coli induced the production of hypoxia-related enzymes (Crespo et al.2016; Jo et al. 2017; Létoffé et al. 2017). Our results, coupled with these previous studies, suggested that biofilm formation was likely beneficial to the pathogenic bacteria to survival under lower oxygen niches.

Unexpectedly, we found that the pathogenic V. parahaemolyticus ATCC17802 was capable of maintaining more rod-shaped cells with no surface shrinkage than the other tested isolates, demonstrating that V. parahaemolyticus ATCC17802 was the most tolerant to the hypoxia. This may attribute to its survial in the host gastrointestinal tract with low oxygen level.

Draft genomes of the four V. parahaemolyticus isolates from aquatic animals were determined (4.91–5.35 Mb), with no plasmid sequenced. Clean single peaks with a typical Poisson distribution in the frequency of observed unique 17-mers within the sequencing data indicated less repetitive DNA in the genomes, consistent with the fewer ISs identified in the V. parahaemolyticus genomes. Among the 4,458–4,924 predicted genes, 735 to 1,046 encoded unknown proteins. A large number of unknown genes have also been reported in seven V. parahaemolyticus genomes of aquatic animal origins in the previous report (Xu et al. 2022a).

V. parahaemolyticus B8-26, B1-21, N2-5, and L7-40 genomes also carried MGEs. Only one identified IS (V. parahaemolyticus L7-40 genome in IS001) was present at the end of the scaffold 84 and scaffold 86 (Table S14), suggesting that the assembled genomes did not result in the loss of identified MGEs (except for this IS001). Of note, the V. parahaemolyticus B1-21 genome had the maximum numbers of INs (n = 11) and prophage gene clusters (n = 2), while the V. parahaemolyticus N2-5, and B8-26 genomes carried the maximum numbers of GIs (n = 14), and ISs (n = 4), respectively. The MGEs, such as GIs (n = 5–9), prophage gene clusters (n = 0–2), INs (n = 1–11), and ISs (n = 0–3) have also been reported in V. parahaemolyticus isolates from six species of aquatic animals (Xu et al. 2022a). Several MGEs were found in V. parahaemolyticus ATCC17802 (Yang et al. 2015) as well, including GIs (n = 2), and INs (n = 7). These results in this study, coupled with the previous reports (Xu et al. 2022a; Yu et al. 2022; Yang et al. 2015), indicated possible HGT via the MGEs during the V. parahaemolyticus evolution. The MGEs, with diverse gene functions identified in this study, can promote ecological niche adaptation, alter gene availability within microbial communities, and drive the bacterial evolution (Tuttle et al. 2022). CRISPR-Cas systems provide adaptive immunity for bacteria to resist foreign DNA invasion (Zaayman et al. 2022). In this study, CRISPR-Cas gene clusters (75–189 bp) were found in the V. parahaemolyticus B1-21 (n = 1), N2-5 (n = 1), L7-40 (n = 3), and ATCC17802 (n = 4) genomes, but none of which carried Cas protein-coding gene, suggesting partial or inactive CRISPR-Cas systems in these V. parahaemolyticus isolates (Fig. S8).

Most types of Vibrio species have an extensive reservoir of virulence genes (Abdelaziz Gobarah et al. 2022). In this study, although no toxic tdh and trh genes (Ceccarelli et al. 2013) were found, virulence-related genes (42–45) were identified in the four V. parahaemolyticus genomes of aquatic animal origins, e.g., ilpA, MAM7, exsA, flagella, Al-2, VpadF, and T3SS1-related genes. T3SSs are important determinants in V. parahaemolyticus pathogenicity (Raval et al. 2021). In our previous study, we also found a number of virulence-related genes (e.g., impI, vgrG, fliCDEFG, flrBC, yscOPQRST, and T6SS-related genes) in V. parahaemolyticus strains recovered from six species of aquatic animals (Xu et al. 2022a). These results suggested possible health risks in consuming these aquatic animals.

In this study, a few antibiotic resistance-related genes, which are responsible for the fluoroquinolone (crp), TET (Tet (35)), and beta-lactam (bla_CARB-18) resistance, were identified in the V. parahaemolyticus genomes, consistent with the observed phenotypes. Moreover, several heavy metal tolerance-related genes (e.g., cusARS, dsbABC, smtA, zntA) were found in these V. parahaemolyticus isolates as well. For example, the smtA and zntA genes, which are involved in bacterial tolerance to Zn, as well as Hg and Zn, respectively (Chen et al. 2021; Lee et al. 2001; Pal et al. 2014), were present in the V. parahaemolyticus B8-26, B1-21, L7-40, and ATCC17802 genomes. V. parahaemolyticus B8-26 also carried the cusARS and dsbABC gene clusters, which are involved in bacterial tolerance to Cu, as well as Cd, Zn, Hg and Cu, respectively. The antibiotic resistance-related genes (such as acrB, catB, hns, qnr, soxR), and heavy metal tolerance-related genes (e.g., copA, cusA, cusR, cusS, zntA) have also been reported in the V. parahaemolyticus genomes derived from six species of aquatic animals (Xu et al. 2022a).

The transcriptomic analyses revealed distinct response strategies to hypoxia by Vibrio parahaemolyticus isolates of clinical and aquatic animal origins. Different degrees of impact on the growth of V. parahaemolyticus isolates were observed, which were triggered by the hypoxic condition. For instance, all DEGs were significantly down-regulated in the PG biosynthesis, as well as ubiquinone and other terpenoid-quinone biosynthesis in V. parahaemolyticus B8-26 (p < 0.05). In response to the hypoxic condition, V. parahaemolyticus L7-40 inhibited the glyoxylate and dicarboxylate metabolism, but up-regulated the glycolysis/gluconeogenesis, and changed the nitrogen metabolism, and the arginine biosynthesis (p < 0.05). The transcriptome data, coupled with the phenotype results in this study, demonstrated that V. parahaemolyticus ATCC17802 was the most tolerant to the hypoxia among the isolates. This isolate down-regulated and/or shut down all the ten changed metabolic pathways to reduce cell viability and maintain cell structure under the hypoxic condition. In particular, the pyrimidine metabolism, fatty acid biosynthesis, LPS biosynthesis, purine metabolism, as well as arginine and proline metabolism were all repressed in V. parahaemolyticus ATCC17802. Fatty acids are essential for biosynthesis of phospholipids and also provide ATP, signaling molecules, and NADPH through β-oxidation (Chen et al. 2022). For example, the DEG encoding a β-ketoacyl-ACP synthase (Vp_17802_11840) in the fatty acid biosynthesis was significantly down-regulated (0.200–fold) in V. parahaemolyticus ATCC17802 (p < 0.05), which play a vital role in cell wall synthesis, biofilm formation and also pathogenesis (Hu et al. 2020). These results provided a novel mechanism for the persistence of the toxic strain in the environment and in the host.

Additionally, some metabolic pathways were altered in the all V. parahaemolyticus isolates under the hypoxic condition, e.g., the inhibited ribosome in V. parahaemolyticus B1-21, N2-5, and ATCC17802 isolates; the repressed aminoacyl-tRNA biosynthesis in V. parahaemolyticus B8-26, and ATCC17802 isolates; the inhibited oxidative phosphorylation in V. parahaemolyticus N2-5, L7-40, and ATCC17802 isolates; and the suppressed citrate cycle in V. parahaemolyticus B8-26, N2-5, and L7-40 isolates (p < 0.05). For instance, the DEGs encoding the fumarate reductase subunits FrdCD, which is involved in maintaining redox balance through regeneration of reduced cofactors during oxygen deficiency conditions (Sędzielewska et al. 2012), were significantly repressed in V. parahaemolyticus B8-26 (Vp_B8-26_15850, 0.107–fold; Vp_B8-26_15855, 0.087–fold), V. parahaemolyticus B1-21 (Vp_B1-21_15850, 0.406–fold; Vp_B1-21_15855, 0.341–fold), V. parahaemolyticus N2-5 (Vp_N2-5_15850, 0.149–fold; Vp_N2-5_15855, 0.230–fold), V. parahaemolyticus ATCC17802 (Vp_17802_15850, 0.248–fold; Vp_17802_15855, 0.113–fold), and V. parahaemolyticus L7-40 (FrdD, Vp_L7-40_15855, 0.417–fold). Most of the DEGs in the oxidative phosphorylation were also significantly repressed in V. parahaemolyticus N2-5, L7-40, and ATCC17802 isolates, which possibly provided indirect evidence for the hypoxic condition maintained in our culture system. The oxidative phosphorylation is the last metabolic pathway of cellular respiration, for the production of ATP, the "energy currency" to sustain life (Wilson 2017).

Taken, V. parahaemolyticus developed multiple molecular strategies to balance the redox state of the cell and the energy conversion, to inhibit the cell wall synthesis, thereby the cell proliferation, to retard the cell growth, thus prolonging the growth cycle and maintaining the cell structure integrated, to enhance the expression of stress-related proteins, structurally stabilizing factors, and efflux RND transporters, in order to efficiently mitigate the cell damage and/or cytotoxicity caused by the hypoxia.

Conclusions

In this study, we investigated, for the first time, the growth of V. parahaemolyticus B8-26, B1-21, N2-5, L7-40, and ATCC17802 isolates under the hypoxic condition. The former four were isolated from S. stritus, C. aurea, O. oratoria, and C. idellus, respectively, while V. parahaemolyticus ATCC17802 was of human clinical origin. These isolates harbored multiple antibiotic and heavy metal resistance phenotypes, and grew optimally in the TSB (3% NaCl, pH 8.5) at 37 °C.

V. parahaemolyticus B8-26, B1-21, N2-5, L7-40, and ATCC17802 isolates had an oxygen concentration-dependent growth mode, and the 10% O₂ condition strongly inhibited the growth of the isolates, when incubated in TSB medium at 37 °C. Unexpectedly, in marked contrast to the normal 21% O₂ condition, the V. parahaemolyticus isolates enhanced biofilm formation under the hypoxic condition for 24 h.

Draft genome sequences of V. parahaemolyticus B8-26 (serotype: O1:K34), B1-21 (serotype: O5:K23), N2-5 (serotype: O4/O11:K4), and L7-40 (serotype: O9: KUT) were determined (4,913,675–5,353,490 bp). Comparative genomic analysis revealed unknown function genes (19.32–23.13%), and MGEs (n = 12–29) in the V. parahaemolyticus genomes. Moreover, V. parahaemolyticus B1-21, and N2-5 had unique genome features, classified into single branches in the phylogenetic tree.

Genome-wide gene expression changes triggered by the hypoxic condition were further examined. Comparative transcriptomic analysis unveiled multiple significantly changed metabolic pathways in the isolates under the 10% O₂ condition for 24 h. V. parahaemolyticus developed multiple molecular strategies to efficiently mitigate the cell damage and/or cytotoxicity caused by the hypoxia.

In addition, V. parahaemolyticus ATCC17802 (serotype: O1:K1) of the clinical origin was the most tolerant to the hypoxia among the test isolates, and down-regulated and/or shut down ten metabolic pathways to reduce cell viability and prolong growth under the stress.

Overall, the results of this study fill prior gaps in V. parahaemolyticus response to the hypoxic condition, and provide a novel mechanism for the toxic V. parahaemolyticus strain to persist in the environment and in the host.