Genomic and Transcriptomic Analysis Reveal Multiple Strategies for the Cadmium Tolerance in Vibrio parahaemolyticus N10-18 Isolated from Aquatic Animal Ostrea gigas Thunberg

The waterborne Vibrio parahaemolyticus can cause acute gastroenteritis, wound infection, and septicemia in humans. Pollution of heavy metals in aquatic environments is proposed to link high incidence of the multidrug-resistant (MDR) pathogen. Nevertheless, the genome evolution and heavy metal tolerance mechanism of V. parahaemolyticus in aquatic animals remain to be largely unveiled. Here, we overcome the limitation by characterizing an MDR V. parahaemolyticus N10-18 isolate with high cadmium (Cd) tolerance using genomic and transcriptomic techniques. The draft genome sequence (4,910,080 bp) of V. parahaemolyticus N10-18 recovered from Ostrea gigas Thunberg was determined, and 722 of 4653 predicted genes had unknown function. Comparative genomic analysis revealed mobile genetic elements (n = 11) and heavy metal and antibiotic-resistance genes (n = 38 and 7). The bacterium significantly changed cell membrane structure to resist the Cd2+ (50 μg/mL) stress (p < 0.05). Comparative transcriptomic analysis revealed seven significantly altered metabolic pathways elicited by the stress. The zinc/Cd/mercury/lead transportation and efflux and the zinc ATP-binding cassette (ABC) transportation were greatly enhanced; metal and iron ABC transportation and thiamine metabolism were also up-regulated; conversely, propanoate metabolism and ribose and maltose ABC transportation were inhibited (p < 0.05). The results of this study demonstrate multiple strategies for the Cd tolerance in V. parahaemolyticus.


Introduction
Vibrio parahaemolyticus is a Gram-negative bacterium that thrives in marine, riverine, and aquaculture environments worldwide [1,2]. The bacterium can cause acute gastroenteritis, wound infection, and septicemia in humans [2]. V. parahaemolyticus was first identified as a foodborne pathogen in Japan in the 1950s [3]. Since then, pathogenic V. parahaemolyticus has been reported in Asian countries and subsequently in Africa, America, and Europe, arguing a pandemic of V. parahaemolyticus worldwide [4]. It was estimated that V. parahaemolyticus is responsible for roughly 35,000 human infection cases each year in the United States [5]. The bacterium has been identified as the leading cause of the foodborne diarrhea disease in China since the 1990s [6]. The crucial virulence determinants in pathogenic V. parahaemolyticus are thermostable-direct hemolysin (TDH) and TDH-related hemolysin (TRH) [7].

Phylogenetic Tree Analysis
Complete gene sequences of 64 V. parahaemolyticus strains were downloaded from the GenBank database (Table S2). Amino acid data sets of single-copy orthologs in V. parahaemolyticus genomes were analyzed using the software OrthoFinder (version 2.2.6) [42]. The FastTree (version 2.1.11) software was used to build a phylogenetic tree using the method and parameters described in our recent research [21].

Determination of Minimum Inhibitory Concentrations (MICs) of Antibiotics and Heavy Metals
The MICs of antibiotics and heavy metals against V. parahaemolyticus N10-18 were measured using the broth dilution testing (microdilution) according to the guidelines of the Clinical and Laboratory Standards Institute (CLSI, M2-A9, 2006), including the CdCl 2 , ZnCl 2, AMP, kanamycin (KAN), and STR (Sinopharm Chemical Reagent Co., Ltd., Shanghai, China). Escherichia coli K-12 was used as a quality control strain in the tests [23].

Stress Conditions
The fresh cell culture of V. parahaemolyticus N10-18 was individually inoculated into the TSB medium supplemented with different concentrations (0, 50,100,200, and 400 µg/mL) of CdCl 2 and then incubated at 37 • C for 48 h. Bacterial growth curves were measured [24]. Bacterial survival rates were calculated using the standard colony counting method [43].

Cell Membrane Permeability, Fluidity, and Surface Hydrophobicity Assays
V. parahaemolyticus N10-18 was incubated in the TSB medium (3% NaCl, pH 8.5) to the mid-LGP at 37 • C. A final concentration of CdCl 2 (50 µg/mL) was added and then incubated at 37 • C for 2 h. The outer cell membrane permeability was examined using the method described by Harman et al. [44]. The N-phenyl-1-naphthylamine (NPN) was purchased from the Shanghai Labtop Bio-Technology Co., Ltd., Shanghai, China. The inner membrane permeability was examined using the method described by Ibrahim et al. [45]. The O-nitrophenyl-β-D galactopyranoside (ONPG) was purchased from Beijing Solarbio Science & Technology Co., Ltd., Beijing, China.
The membrane fluidity assay was performed using the method described by Voss and Montville [46]. The cell surface hydrophobicity assay was performed using the method described by Yan et al. [47]. The n-hexadecane was purchased from China National Pharmaceutical Group Corporation Co., Ltd., Shanghai, China.

Scanning Electron Microscope (SEM) Assay
The SEM assay was performed according to the method described previously [48]. Briefly, a final concentration (50 µg/mL) of CdCl 2 was added into V. parahaemolyticus N10-18 culture grown in the TSB medium (pH 8.5, 3% NaCl) at mid-LGP and then continuously incubated at 37 • C for 2 h. An amount of 1.5 mL of the cell suspension was collected, washed, dehydrated, dried, and gold-covered by cathodic spraying and observed using the thermal field emission SEM (Hitachi, SU5000, Tokyo, Japan) with accelerating voltages of 5-10 kV [48].

Illumina RNA Sequencing and Analysis
V. parahaemolyticus N10-18 was incubated in the TSB medium (pH 8.5, 3% NaCl) to the mid-LGP at 37 • C. A final concentration of CdCl 2 (50 µg/mL) was added and then incubated at 37 • C for 2 h. Controls were cultures also exposed to no cadmium for the same time period and collected as treatments. The bacterial cells were collected by centrifugation and subjected for the Illumina RNA sequencing. The RNA extraction and quality control, sequencing library construction, and Illumina sequencing were conducted by Shanghai Majorbio Bio-pharm Technology Co., Ltd., Shanghai, China, using Illumina HiSeq 2500 platform (Illumina, Santiago, CA, USA). Three replicates were conducted for each sample.
Expression of each gene was calculated, and differentially expressed genes (DEGs) were defined and used for gene set enrichment analysis (GSEA) as described previously [25]. Representative DEGs were examined using real-time reverse-transcription PCR (RT-qPCR) assay [25,48].

Statistical Analysis
The data were analyzed using the SPSS software (version 22, IBM, Armonk, NY, USA). All tests in this study were conducted in triplicate.

Results
3.1. Genotype and Phenotype of V. parahaemolyticus N10-18 V. parahaemolyticus N10-18 isolate was recovered from O. gigas Thunberg [23]. The bacterium tested negative for the toxic tdh and trh genes but positive for the species-specific gene tlh [49]. The results also showed that V. parahaemolyticus N10-18 was tolerant to the heavy metals Cd 2+ and Zn 2+ , as well as the antimicrobial agents AMP, KAN, and STR (Table S1).

Genome Features of V. parahaemolyticus N10-18
The ANI value of the V. parahaemolyticus N10-18 genome was determined, which was higher (98.22%) than the threshold (94-96%) for species determination [50]. The draft genome sequence of V. parahaemolyticus N10-18 was determined using the Illumina Hiseq × 10 sequencing technique (Figure 1), and approximately 438,181 clean single reads were obtained. The assembly generated 70 scaffolds with a sequencing depth (on average) of 319.27-fold. V. parahaemolyticus N10-18 showed a clean single peak in the frequency of observed unique 17-mers within the sequencing data and varied as a typical Poisson distribution, suggesting less repetitive DNA in the V. parahaemolyticus N10-18 genome ( Figure S1). The obtained V. parahaemolyticus N10-18 genome sequence has been deposited in the GenBank database under the assigned accession number JALGSE000000000. The V. parahaemolyticus N10-18 genome contained transposase genes (n = 10) and MGEs, including GIs (n = 2), INs (n = 8), and ISs (n = 1), suggesting possible horizontal gene transfer (HGT) mediated by the MGEs during the V. parahaemolyticus N10-18 genome evolution. The identified MGEs were absent from the ends of the scaffolds (Table S3), which indicated that the draft genome contained all such elements. Circles from the inside to outside: GC contents (outward part means higher than average, inward part means lower than average); GC skew (purple value is greater than zero, green value is less than zero); the reference genome of V. parahaemolyticus RIMD2210633 (GenBank accession numbers: NC_004603.1 and NC_004605.1) and V. parahaemolyticus N10-18 genome (GenBank accession no. JALGSE000000000), respectively; and CDSs on the positive and negative chains, respectively.

Serotype and ST of V. parahaemolyticus N10-18
The BLAST analysis of the antigen gene loci revealed that the V. parahaemolyticus N10-18 genome contained the O antigen loci orf16/wvdB and specific loci wzc for K4 polymorphic sites [51], indicating that the serotype of V. parahaemolyticus N10-18 was O4/O11:K4. Additionally, the ST by the MLST analysis showed that the bacterium belonged to the ST-499.

Phylogenetic Relatedness of V. parahaemolyticus N10-18
Approximately 1485 homologous single-copy amino acid sequences were identified from 64 V. parahaemolyticus genomes available in the GenBank database and the V. para- Figure 1. Genome circle maps of V. parahaemolyticus N10-18. (A,B) represent the larger and smaller chromosomes of V. parahaemolyticus N10-18, respectively. Circles from the inside to outside: GC contents (outward part means higher than average, inward part means lower than average); GC skew (purple value is greater than zero, green value is less than zero); the reference genome of V. parahaemolyticus RIMD2210633 (GenBank accession numbers: NC_004603.1 and NC_004605.1) and V. parahaemolyticus N10-18 genome (GenBank accession no. JALGSE000000000), respectively; and CDSs on the positive and negative chains, respectively.
The obtained genome size of V. parahaemolyticus N10-18 was 4,910,080 bp with 45.46% of the GC content. A total of 4653 genes were predicted, among which approximately 4565 coded for proteins. Remarkably, approximately 722 proteins-coding genes had unknown  (Table 1). The V. parahaemolyticus N10-18 genome contained transposase genes (n = 10) and MGEs, including GIs (n = 2), INs (n = 8), and ISs (n = 1), suggesting possible horizontal gene transfer (HGT) mediated by the MGEs during the V. parahaemolyticus N10-18 genome evolution. The identified MGEs were absent from the ends of the scaffolds (Table S3), which indicated that the draft genome contained all such elements.

Serotype and ST of V. parahaemolyticus N10-18
The BLAST analysis of the antigen gene loci revealed that the V. parahaemolyticus N10-18 genome contained the O antigen loci orf16/wvdB and specific loci wzc for K4 polymorphic sites [51], indicating that the serotype of V. parahaemolyticus N10-18 was O4/O11:K4. Additionally, the ST by the MLST analysis showed that the bacterium belonged to the ST-499.

INs
Mobile INs are prevalent in human-dominated ecosystems with prolonged exposure to selective agents such as detergents, antibiotics, and heavy metals [53]. INs are generally classified according to integrase genes intI 1, intI 2, intI 3, and intI 4 into type I, type II, type III, and super integron, respectively [54].  reading frame (orf ) and a recombination site (attC) necessary for integration. They can exist free as circular molecules or mobilized in INs [55].

INs
Mobile INs are prevalent in human-dominated ecosystems with prolonged exposure to selective agents such as detergents, antibiotics, and heavy metals [53]. INs are generally classified according to integrase genes intI 1, intI 2, intI 3, and intI 4 into type I, type II, type III, and super integron, respectively [54]. In this study, eight INs (IN 1 to IN 8) were identified in the V. parahaemolyticus N10-18 genome, which ranged from 910 bp to 227,599 bp and carried 2 to 210 genes. Of these, there was one complete IN (IN 1) and seven gene cassettes (IN 2 to IN 8) ( Figure 4). Typically, gene cassettes consist of a promoterless open reading frame (orf) and a recombination site (attC) necessary for integration. They can exist free as circular molecules or mobilized in INs [55].
In this study, the complete IN 1 (2566 bp) contained a hypothetical protein-encoding gene (Vp_N10_18_2516) and an integrase gene IntI 4 (Vp_N10_18_2515). The latter showed sequence identity (99.38%) with the super IN IntI 4 (NR reference sequence: AHI99301.1) [56], which indicated that IN

ISs
A single short IS can transfer one or more resistance-related genes in Gram-negative bacteria and affect bacterial resistance phenotype [21,22]. In this study, only one IS110 (1327 bp) was identified in the V. parahaemolyticus N10-18 genome, encoding a IS110 family transposase (Table S3).

Heavy Metal and Antibiotic Resistance-Associated Genes
Approximately 38 heavy metal tolerance-associated genes were identified in the V.

ISs
A single short IS can transfer one or more resistance-related genes in Gram-negative bacteria and affect bacterial resistance phenotype [21,22]. In this study, only one IS110 (1327 bp) was identified in the V. parahaemolyticus N10-18 genome, encoding a IS110 family transposase (Table S3).

Heavy Metal and Antibiotic Resistance-Associated Genes
Approximately 38 heavy metal tolerance-associated genes were identified in the V. parahaemolyticus N10-18 genome by the BLAST analysis (Table 2). For example, the cadC gene and dsbABC gene cluster, which are responsible for the bacterial tolerance to Cd, Zn, and Pb, as well as Cd, Zn, Hg, and Cu, respectively [65], were present in the V. parahaemolyticus N10-18 genome. Moreover, the zntAR, znuABC, zur, and smtA genes for the Zn and Hg tolerance [65][66][67] were also identified ( Table 2). These results were consistent with the heavy metal tolerance phenotype of V. parahaemolyticus N10-18. Antimicrobial resistance-associated genes (n = 7) also existed in the V. parahaemolyticus N10-18 genome (Table 2), such as an elongation factor Tu (tuf ) [68], a cAMP-activated global transcriptional regulator CRP (crp) [69], a DNA-directed RNA polymerase subunit beta (rpoB) [70], and a hexose-6-phosphate (uhpT) [71], tet34 and tet35 [72], and β-lactamase (bla CARB-21 ), consistent with the MDR phenotype of V. parahaemolyticus N10-18.  (Table S1). Remarkably, the observed MICs of Cd 2+ and Zn 2+ were 400 µg/mL and 1600 µg/mL, respectively. Given that Zn 2+ is essential for the growth and development of aquatic animals and often used as feed additives, the survival of V. parahaemolyticus N10-18 to resist the high level of Cd 2+ was further investigated in this study.
Growth curves of V. parahaemolyticus N10-18 at different concentrations of CdCl 2 were determined in the TSB medium at 37 • C. As shown in Figure 5 (Table S1). Remarkably, the observed MICs of Cd 2+ and Zn 2+ were 400 μg/mL and 1600 μg/mL, respectively. Given that Zn 2+ is essential for the growth and development of aquatic animals and often used as feed additives, the survival of V. parahaemolyticus N10-18 to resist the high level of Cd 2+ was further investigated in this study.
Growth curves of V. parahaemolyticus N10-18 at different concentrations of CdCl2 were determined in the TSB medium at 37 °C. As shown in Figure 5, at the concentration of 400 μg/mL of Cd 2+ , the growth of V. parahaemolyticus N10-18 was fully inhibited. At 200 μg/mL and 100 μg/mL of Cd 2+ , the bacterial growth was retarded, showing a longer lag phase of 32 h and 6 h, respectively. Moreover, the bacterial biomass reached the maximum with OD600 nm values of 1.178 and 1.216 at 48 h and 28 h, respectively. At 50 μg/mL of Cd 2+ , only a slight decrease in growth was observed, when compared with the control (0 μg/mL of Cd 2+ ). Under this treatment, the observed fatality rate of V. parahaemolyticus N10-18 was 10.73%.  Bacterial cell membrane permeability and fluidity and cell surface hydrophobicity are key parameters of the cell membrane responding to environmental challenges [80,81]. Cd is a heavy metal whose cations often cause toxicity to both eukaryotic and prokaryotic cells even at low concentrations [82]. In this study, the outer cell membrane permeability was examined using the NPN probe. As shown in Figure 6A, the probe fluorescence intensity showed an overall downward trend after treatment with 50 µg/mL of Cd 2+ for 4 h. Additionally, we used the ONPG as a probe to examine the inner cell membrane permeability, and no significant difference in the inner membrane permeability was also observed after the Cd 2+ stress for 1.5 h compared with the control group (p > 0.05). However, the extended treatment time (2.0 to 4.0 h) increased the bacterial inner membrane permeability by 2.04to 4.96-fold (p < 0.05) ( Figure 6B). was examined using the NPN probe. As shown in Figure 6A, the probe fluorescence intensity showed an overall downward trend after treatment with 50 μg/mL of Cd 2+ for 4 h. Additionally, we used the ONPG as a probe to examine the inner cell membrane permeability, and no significant difference in the inner membrane permeability was also observed after the Cd 2+ stress for 1.5 h compared with the control group (p > 0.05). However, the extended treatment time (2.0 to 4.0 h) increased the bacterial inner membrane permeability by 2.04-to 4.96-fold (p < 0.05) ( Figure 6B). Cytoplasmic membrane fluidity also influences the ability of most compounds (nutrients and antibiotics) and ions to cross the bacterial cytoplasmic membrane by diffusion and active transport [80]. As shown in Figure 6C, cell membrane fluidity of V. parahaemolyticus N10-18 was significantly decreased by 1.07-fold after being treated with 50 μg/mL of Cd 2+ for 2 h, compared with the control group. Cell surface hydrophobicity is crucial in bacterial adhesion to abiotic and biological surfaces [83]. As shown in Figure 6D, cell surface hydrophobicity of V. parahaemolyticus N10-18 was significantly increased by 1.47-fold after being treated with 50 μg/mL of Cd 2+ for 2 h (p < 0.05).
Cell structure changes of V. parahaemolyticus N10-18 under the Cd 2+ (50 μg/mL) stress were also observed by the SEM analysis. As shown in Figure 7, the treatment with Cd 2+ (50 μg/mL) for 2 h resulted in the cell surface shrinking of certain V. parahaemolyticus N10-18 cells compared to the control group. Cytoplasmic membrane fluidity also influences the ability of most compounds (nutrients and antibiotics) and ions to cross the bacterial cytoplasmic membrane by diffusion and active transport [80]. As shown in Figure 6C, cell membrane fluidity of V. parahaemolyticus N10-18 was significantly decreased by 1.07-fold after being treated with 50 µg/mL of Cd 2+ for 2 h, compared with the control group. Cell surface hydrophobicity is crucial in bacterial adhesion to abiotic and biological surfaces [83]. As shown in Figure 6D, cell surface hydrophobicity of V. parahaemolyticus N10-18 was significantly increased by 1.47-fold after being treated with 50 µg/mL of Cd 2+ for 2 h (p < 0.05).
Cell structure changes of V. parahaemolyticus N10-18 under the Cd 2+ (50 µg/mL) stress were also observed by the SEM analysis. As shown in Figure 7, the treatment with Cd 2+ (50 µg/mL) for 2 h resulted in the cell surface shrinking of certain V. parahaemolyticus N10-18 cells compared to the control group.

The Major Changed Metabolic Pathways Medicated by the Cd 2+ (50 µg/mL) Stress in V. parahaemolyticus N10-18
Based on the obtained results, V. parahaemolyticus N10-18 grown at the mid-LGP in the TSB medium at 37 • C was treated with the Cd 2+ (50 µg/mL) for 2 h, and gene expression changes at the global genome level of V. parahaemolyticus N10-18 induced by the Cd 2+ stress were determined using the Illumina HiSeq 2500 sequencing technology.

The Major Changed Metabolic Pathways Medicated by the Cd 2+ (50 μg/mL) Stress in V. parahaemolyticus N10-18
Based on the obtained results, V. parahaemolyticus N10-18 grown at the mid-LGP in the TSB medium at 37 °C was treated with the Cd 2+ (50 μg/mL) for 2 h, and gene expression changes at the global genome level of V. parahaemolyticus N10-18 induced by the Cd 2+ stress were determined using the Illumina HiSeq 2500 sequencing technology. Approximately 8.3% (377 of 4565 genes) of the bacterial genes were differentially expressed under the treatment, when compared to the control group. Of these, 217 DEGs showed higher transcriptional levels (fold change ≥ 2.0), whereas 160 were significantly down-regulated (fold change ≤ 0.5) ( Figure 8A). Approximately seven significantly altered metabolic pathways were identified in V. parahaemolyticus N10-18, including the ATPbinding cassette (ABC) transporters, propanoate metabolism, benzoate degradation, thiamine metabolism, fat digestion and absorption, quorum sensing (QS), and pathogenic E. coli infection ( Figure 8B, Table 3).
Remarkably, the expression of approximately 28 DEGs of the ABC transporters was significantly changed at the transcription level (0.061-to 11.609-fold) (p < 0.05). Of these, the DEGs encoding the maltose and ribose transporters and some amino acid transporters were significantly inhibited (0.061-to 0.500-fold). For example, the malEK gene cluster (Vp_N10_18_1557, Vp_N10_18_1556), which encoded a maltose ABC transporter substrate-binding protein MalE and a maltose/maltodextrin import ATP-binding protein MalK, respectively, was significantly down-regulated. The rbsBCD gene cluster (Vp_N10_18_3026, Vp_N10_18_3025, and Vp_N10_18_3023), which encoded a ribose ABC transporter substrate-binding protein RbsB, a ribose ABC transporter permease, and a Dribose pyranase, respectively, was significantly down-regulated as well (p < 0.05). Additionally, the expression of the livH gene (Vp_N10_18_2959), which encoded a branchedchain amino acid ABC transporter permease, was remarkably down-regulated (0.061-fold) at the transcriptional level.   In the propanoate metabolism, expression of approximately seven DEGs was also significantly inhibited (0.069-to 0.438-fold) (p < 0.05), e.g., the prpCEF, and acnD genes (Vp_N10_18_0015, Vp_N10_18_0011, Vp_N10_18_0012, and Vp_N10_18_0013) involved in the 2-methylcitric acid cycle (2-MCC). The 2-MCC in the propionate catabolic pathway is used to oxidize the Cα methylene of propionate to a keto group yielding pyruvate [84]. These results indicated that V. parahaemolyticus N10-18 greatly reduced the branchedchain amino acid transportation, inhibited the maltose and ribose transportation, and inactivated to utilize the propionic acid as a carbon source under the Cd 2+ (50 μg/mL) stress.
Conversely, the DEGs encoding the Zn/Cd/Hg/Pb-transporting ATPase (zntA, Vp_N10_18_0526) and heavy metal efflux resistance-nodulation-cell division (RND) trans-  Remarkably, the expression of approximately 28 DEGs of the ABC transporters was significantly changed at the transcription level (0.061-to 11.609-fold) (p < 0.05). Of these, the DEGs encoding the maltose and ribose transporters and some amino acid transporters were significantly inhibited (0.061-to 0.500-fold). For example, the malEK gene cluster (Vp_N10_18_1557, Vp_N10_18_1556), which encoded a maltose ABC transporter substrate-binding protein MalE and a maltose/maltodextrin import ATP-binding protein MalK, respectively, was significantly down-regulated. The rbsBCD gene cluster (Vp_N10_18_3026, Vp_N10_18_3025, and Vp_N10_18_3023), which encoded a ribose ABC transporter substratebinding protein RbsB, a ribose ABC transporter permease, and a D-ribose pyranase, respectively, was significantly down-regulated as well (p < 0.05). Additionally, the expression of the livH gene (Vp_N10_18_2959), which encoded a branched-chain amino acid ABC transporter permease, was remarkably down-regulated (0.061-fold) at the transcriptional level.
In the propanoate metabolism, expression of approximately seven DEGs was also significantly inhibited (0.069-to 0.438-fold) (p < 0.05), e.g., the prpCEF, and acnD genes (Vp_N10_18_0015, Vp_N10_18_0011, Vp_N10_18_0012, and Vp_N10_18_0013) involved in the 2-methylcitric acid cycle (2-MCC). The 2-MCC in the propionate catabolic pathway is used to oxidize the Cα methylene of propionate to a keto group yielding pyruvate [84]. These results indicated that V. parahaemolyticus N10-18 greatly reduced the branched-chain amino acid transportation, inhibited the maltose and ribose transportation, and inactivated to utilize the propionic acid as a carbon source under the Cd 2+ (50 µg/mL) stress.
Interestingly, the DEGs encoding the T3SS needle filament protein VscF (Vp_N10_18 _0060) and glyceraldehyde-3-phosphate dehydrogenase (gapA, Vp_N10_18_3876) in pathogenic E. coli infection were significantly up-regulated by 5.836-and 2.086-fold, respectively (p < 0.05). The VscF in T3SS1 helps to translocate VPA0226 that can be secreted into the host cell cytoplasm via T3SS1 in V. parahaemolyticus [58]. The gapA gene was only expressed under certain stress conditions, and overproduction of GapA led to increased resistance to H 2 O 2 in Lactococcus lactis MG1363 [87]. Additionally, expression of representative DEGs was examined by the RT-PCR assay, and the resulting data were generally consistent with the transcriptomic analysis (Tables S6 and S7).
Taken together, under the Cd 2+ (50 µg/mL) stress, V. parahaemolyticus N10-18 employed multiple strategies for efficient transportation and exocytosis of Cd 2+ to alleviate its cytotoxicity: (1) greatly enhanced the Zn/Cd/Hg/Pb-transportation and efflux; (2) upregulated metal and iron ABC transportation; (3) enhanced the biosynthesis of TPP in the thiamine metabolism; (4) up-regulated the expression of stress-related proteins, such as GapA, and structurally stabilizing factor arginine; (5) conversely, greatly reduced the branched-chain amino acid transportation; (6) inhibited the maltose and ribose ABC transportation; and (7) down-regulated the propanoate metabolism, in order to survive in the adverse niche.

Discussion
The pollution of heavy metals in aquatic environments has led to heavy metal residues in aquatic products, which poses a huge hidden danger in human health [88,89]. Cd is classified into Group 1 as carcinogenic to humans by the IARC [14]. This toxic element possibly results in short-term or long-term disorders in the body, such as degenerative bone disease, kidney dysfunction, lung injuries, disorders in Zn and Cu metabolism, and cancer [14]. The heavy metal pollution has also been supposed to link high incidence of the MDR V. parahaemolyticus, which is a challenging issue in the clinical treatment. This study was the first to characterize the MDR V. parahaemolyticus N10-18 with high tolerance to Cd 2+ and Zn 2+ . The observed MIC values of Cd 2+ , Zn 2+ , AMP, KAN, and STR against V. parahaemolyticus N10-18 were 400 µg/mL, 1600 µg/mL, 50,000 µg/mL, 128 µg/mL, and 128 µg/mL, respectively, which suggested possible antibiotic and heavy metal pollution in the aquaculture environment of O. gigas Thunberg, consistent with previous reports [8,15,21,23].
In this study, the draft genome sequence (4,910,080 bp) of V. parahaemolyticus N10-18 was determined using the Illumina Hiseq × 10 sequencing technique. Some repeats were located at the end of scaffolds (n = 23, <1.1 Kb) (Table S5), indicating that the genome assembly was incomplete and contained some gaps, due to limitations of the second-generation Illumina short-read sequencing technique. Approximately 4653 genes were predicted, of which 722 encoded unknown proteins. In our previous study, we found a number of unknown-function genes in V. parahaemolyticus isolates from the six species of aquatic animals (P. undulate, P.a viridis, M.veneriformis, A. nobilis, C. auratu, and L. vannamei) [21]. These results demonstrated genome variation and plasticity of V. parahaemolyticus in aquatic animals.
Intercellular transmissibility of MGEs may have constituted important driving forces in V. parahaemolyticus genome evolution and formation of ecotypes and speciation [90]. In this study, we identified 11 MGEs in the V. parahaemolyticus N10-18 genome. It cannot exclude that additional MGEs may be located in the gaps between scaffolds [21]. In this study, we found two GIs that contained 28 genes in the V. parahaemolyticus N10-18 genome, which facilitated the bacterium to better fit into the niche. For example, GI 1 contained the genes encoding cold-shock proteins (CSPs) (Vp_N10_18_3249, Vp_N10_18_3253). The CSPs served as global gene expression regulators to respond to different stress conditions [22].
INs allow bacteria to capture, stockpile, express, and exchange genes embedded within gene cassettes [91]. In this study, one super IN and seven incomplete INs were identified in the V. parahaemolyticus N10-18 genome. Although approximately 33.20% of the INs-carrying genes (n = 85) encoded unknown proteins, the identified INs endowed the bacterium with diverse environmental adaptability. For instance, several gene cassettes were found to carry virulence-related genes, such as the GNAT super family proteins, prevent-host-death family proteins, and plasmid stabilization proteins. Of these, the GNAT (Vp_N10_18_4551) belongs to type II toxin of toxin-antitoxin systems. The GNAT toxin blocks protein translation by acetylating the amino group of charged tRNAs, thus preventing tRNA from participating in peptidyl ribosomal transferase [92].
Some V. parahaemolyticus isolates lacking the tdh and/or trh genes are also highly cytotoxic to human gastrointestinal cells, which indicates that other virulence factors exist. In this study, we found potential virulence-related genes (n = 45) in the V. parahaemolyticus N10-18 genome, e.g., ilpA, MAM7, exsACD, gmhA, gmd, kdsA, and T3SS1-related genes, which are involved in bacterial adhesion or epithelial cell invasion. For example, the ilpA gene encodes an immunogenic lipoprotein A in Vibrio vulnificus, an adhesion molecule that can induce cytokine production in human immune cells [60]. MAM7 enables Gramnegative pathogens to establish high-affinity binding to host cells during the early stages of infection and is crucial for the delivery of virulence factors into hosts [61]. In this study, we identified 36 genes in T3SS1 in the V. parahaemolyticus N10-18 genome, which are important determinants of the pathogenicity of V. parahaemolyticus. Of these genes, VscCD genes not only activated bacterial resistance to acid stress, H 2 O 2 , and antibiotics but also enhanced the colonization ability and pathogenicity of Vibrio harveyi [93]. These results suggested a health risk in consuming O. gigas Thunberg contaminated by V. parahaemolyticus N10-18.
Bacterial MDR is regarded as an emerging pollutant in different food production avenues including aquaculture [94]. Resistance factors have been reported in pathogenic bacteria [68][69][70][71][72][95][96][97][98]. In this study, we identified seven antibiotic-resistance-related genes in the V. parahaemolyticus N10-18 genome, consistent with the observed MDR phenotype of the bacterium. Remarkably, 38 heavy metal tolerance-associated genes existed in V. parahaemolyticus N10-18. For instance, the gene (Vp_N10_18_3808, GI 2) encoding a short-chain dehydrogenase/reductase SDR family member was identified, which functions in the Cd 2+ stress in Pleurotus eryngii [99]. Moreover, the cadC, dsbABC, zntAR, znuABC, zur, and smtA genes related to the Cd and Zn resistance were also present in the V. parahaemolyticus N10-18 genome, consistent with the high Cd and Zn tolerance phenotype of the bacterium. For instance, the dsbABC gene cluster involved in the degradation of pyrimidine ribonucleosides was found to be related to the resistance and absorbing of Cd in Enterococcus faecalis LZ-11, which was isolated from Lanzhou reach of the Yellow River in China [100]. The diversity of resistance genes, gene variance, and selective pressure from the environment may result in the difference between resistance phenotype and resistance genotype.
In this study, the constructed phylogenetic tree showed that the 65 V. parahaemolyticus genomes were clustered into four large clusters, among which V. parahaemolyticus N10-18 fell into a single sub-branch in Group 4b. The bacterium is located phylogenetically distant from the other V. parahaemolyticus strains originating in aquatic animals but showed the closest evolutionary distance with 8 V. parahaemolyticus strains isolated from Homo sapiens between 1998 and 2015 in the USA. Until 1996, V. parahaemolyticus infection cases were sporadic, occurred in certain countries, and could be related to diverse serovars [101]. Location and isolation time of V. parahaemolyticus strains were not associated with evolutionary taxa, suggesting that the widespread global trade in aquatic products over the past 30 years may have contributed to the cross-regional spread of the pathogen, leading to an increased risk of edible aquatic animals.
Based on the findings in this study, the molecular mechanism underlying the heavy metal Cd 2+ tolerance of V. parahaemolyticus N10-18 was further explored. Under the Cd 2+ (50 µg/mL) stress, the bacterium significantly changed cell membrane permeability and fluidity and cell surface hydrophobicity (p < 0.05). Cell osmotic changes have been disclosed as stressors that can affect biophysical properties and the composition of the membrane and consequently transport mechanisms (permeability) and cell shape and integrity [80]. In this study, after the Cd 2+ (50 µg/mL) treatment, the cell surface of V. parahaemolyticus N10-18 was observed shrinking to a certain extent.
Comparative transcriptomic analysis revealed seven significantly altered metabolic pathways in V. parahaemolyticus N10-18 under the Cd 2+ stress. Remarkably, the DEGs encoding the Zn/Cd/Hg/Pb-transporting ATPase (zntA), and heavy metal efflux RND transporter of the CusA/CzcA family (cusA) were greatly up-regulated by 23.639-and 8.649fold, respectively (p < 0.05). The zntA gene was originally described as a Zn-transporting ATPase in E. coli, but it also confers resistance to Cd [66]. RND efflux pumps are essential for the expulsion of a plethora of potentially small lethal agents or compounds such as detergents, solvents, heavy metals, antibiotics, and toxic secondary metabolites [102]. The CusC(F)BA complex exports copper (I) and silver (I) and mediates resistance to these two metal ions in E. coli [102]. Interestingly, in this study, the ABC transporter encoded by the znuABC genes for high-affinity Zn 2+ uptake [66,67] was also highly increased at the transcriptional level (2.594-to 11.609-fold) (p < 0.05). This was consistent with the high tolerance of V. parahaemolyticus N10-18 to Zn 2+ . Given that Cd is chemically similar to Zn, both of which belong to the IIB transition elements [67], our result provided evidence for the expulsion of Cd 2+ via the Zn 2+ channels in V. parahaemolyticus N10-18. On the other hand, the DEGs (afuA and fhuB) encoding the ion and metal transporters were also up-regulated (p < 0.05). For instance, the fhuB gene encoded a Fe 3+ -hydroxamate ABC transporter permease FhuB. Iron (III) hydroxamate transport across the cytoplasmic membrane is catalyzed by the very hydrophobic FhuB protein and the membrane-associated FhuC protein [103]. These results indicated highly enhanced expulsion of Cd 2+ by V. parahaemolyticus N10-18 to alleviate its cytotoxicity.
In the QS, expression of five DEGs was also significantly up-regulated by 2.14-to 9.727-fold (p < 0.05). Additionally, the pcaC gene (Vp_N10_18_2971) encoding the carboxymuconolactone decarboxylase family protein in benzoate degradation and the atoB gene (Vp_N10_18_2988) encoding the thiolase family protein in the fat digestion and absorption were also significantly enhanced by 2.003-and 4.215-fold, respectively (p < 0.05). These results suggested possibly increased substance absorption for energy conservation and stringent response regulation in V. parahaemolyticus N10-18 under the Cd 2+ (50 µg/mL) stress.
In contrast, the DEGs involved in the branched-chain amino acid transportation and maltose and ribose transportation were significantly repressed (0.061-to 0.500-fold) (p < 0.05). Meanwhile, the propanoate metabolism was also significantly inhibited (0.069to 0.438-fold) (p < 0.05). These results suggested possible repressed energy consumption and nucleotide and ribosome biosynthesis under the Cd 2+ adverse condition. It will be interesting to further investigate the DEGs using proteomic, cell, and animal mode techniques and methods in the future research.

Conclusions
This study was the first to characterize the MDR V. parahaemolyticus N10-18 with high tolerance to Cd 2+ and Zn 2+ (MIC S : 400 µg/mL and 1600 µg/mL) using genomic and transcriptomic techniques. The draft genome sequence (4,910,080 bp) of V. parahaemolyticus N10-18 was determined, and 722 of 4653 predicted genes had unknown function. Comparative genomic analyses revealed MGEs, including GIs (n = 2), INs (n = 8), and ISs (n = 1). Heavy metal and antibiotic-resistance genes (n = 38 and 7) and virulence-associated genes (n = 45) were also found in the V. parahaemolyticus N10-18 genome. The bacterial growth was slightly decreased under the 50 µg/mL of Cd 2+ . V. parahaemolyticus N10-18 significantly changed cell membrane permeability and fluidity and surface hydrophobicity under the Cd 2+ (50 µg/mL) stress (p < 0.05). Meanwhile, comparative transcriptomic analysis revealed seven significantly altered metabolic pathways. Under the Cd 2+ stress, V. parahaemolyticus N10-18 employed multiple strategies for efficient transportation and exocytosis of Cd 2+ to alleviate its cytotoxicity, including greatly enhanced Zn/Cd/Hg/Pb transportation and efflux and significantly up-regulated metal and iron ABC transportation, thiamine metabolism, and stress-related protein expression (e.g., GapA and arginine); in contrast, it greatly reduced the branched-chain amino acid transportation and significantly inhibited the maltose and ribose ABC transportation and propanoate metabolism, in order to resist and survive in the adverse Cd 2+ environment. The results also provided evidence for the expulsion of Cd 2+ via the Zn 2+ channels in V. parahaemolyticus N10-18. Overall, the results of this study enriched genome data of V. parahaemolyticus from aquatic animals and revealed multiple strategies for the cadmium tolerance in the leading seafood-borne pathogen worldwide.