Genomic and Transcriptomic Analyses Reveal Multiple Strategies for Vibrio parahaemolyticus to Tolerate Sub-Lethal Concentrations of Three Antibiotics

Vibrio parahaemolyticus can cause acute gastroenteritis, wound infections, and septicemia in humans. The overuse of antibiotics in aquaculture may lead to a high incidence of the multidrug-resistant (MDR) pathogen. Nevertheless, the genome evolution of V. parahaemolyticus in aquatic animals and the mechanism of its antibiotic tolerance remain to be further deciphered. Here, we investigated the molecular basis of the antibiotic tolerance of V. parahaemolyticus isolates (n = 3) originated from shellfish and crustaceans using comparative genomic and transcriptomic analyses. The genome sequences of the V. parahaemolyticus isolates were determined (5.0–5.3 Mb), and they contained 4709–5610 predicted protein-encoding genes, of which 823–1099 genes were of unknown functions. Comparative genomic analyses revealed a number of mobile genetic elements (MGEs, n = 69), antibiotic resistance-related genes (n = 7–9), and heavy metal tolerance-related genes (n = 2–4). The V. parahaemolyticus isolates were resistant to sub-lethal concentrations (sub-LCs) of ampicillin (AMP, 512 μg/mL), kanamycin (KAN, 64 μg/mL), and streptomycin (STR, 16 μg/mL) (p < 0.05). Comparative transcriptomic analyses revealed that there were significantly altered metabolic pathways elicited by the sub-LCs of the antibiotics (p < 0.05), suggesting the existence of multiple strategies for antibiotic tolerance in V. parahaemolyticus. The results of this study enriched the V. parahaemolyticus genome database and should be useful for controlling the MDR pathogen worldwide.


Introduction
Vibrio parahaemolyticus, a halophilic Gram-negative bacterium, is typically present in marine and estuarine environments [1].Consuming raw or improperly cooked seafood contaminated with pathogenic V. parahaemolyticus can result in acute gastroenteritis, and in severe cases, even death [2].This bacterium was initially recognized as a pathogen in 1950; since then, it has become a major foodborne pathogen worldwide [1,3].Pathogenic V. parahaemolyticus produces thermostable direct hemolysin (TDH) and TDH-associated hemolysin (TRH), which are molecular markers of its pathogenicity [4,5].
V. parahaemolyticus is commonly detected in fish and shellfish worldwide [6][7][8].For example, Vu et al. collected 120 seafood samples from different traditional markets in Hanoi, Vietnam between May and October of 2020 [7].They found that V. parahaemolyticus was present in 58.33% of the samples, with shrimp samples having the highest detection rate at 86.67%, followed by 53.33% of fish samples (n = 30), 53.33% of squid samples (n = 30) and 40% of oyster samples (n = 30) [7].Antibiotics effectively eliminated infectious diseases in aquaculture caused by pathogenic bacteria [9].However, the excessive use of antibiotics in the clinical and aquaculture industries have led to the emergence and evolution of antibioticresistant V. parahaemolyticus over the past few decades [10].V. parahaemolyticus is commonly resistant to ampicillin (AMP), kanamycin (KAN), and streptomycin (STR) [8,11].For example, Lopatek et al. reported that V. parahaemolyticus isolates were found in 595 seafood samples with a detection rate of 17.5% (n = 104).Among these isolates, 75.0% were resistant to AMP, and 68.3% were resistant to STR.Additionally, the majority of the isolates (55.8%) showed resistance to two classes of antimicrobials, primarily AMP and STR [8].Multidrug-resistant (MDR) pathogens in aquaculture environments have the potential to enter aquatic organisms through the food chain and subsequently pose a significant risk to both aquaculture and human health [12].The frequent presence of MDR V. parahaemolyticus has become a critical public health concern [13].
Previous studies have shown that mobile genetic elements (MGEs) facilitate the accumulation and spread of antimicrobial resistance genes through horizontal gene transfer (HGT) [14].With the continuous breakthrough of genome sequencing technology [15], 1740 V. parahaemolyticus isolates have been sequenced so far (GenBank database, https:// www.ncbi.nlm.nih.gov/;accession date: 29 January 2022), of which the complete genomes of 64 V. parahaemolyticus isolates are available in the GenBank database.This provided the possibility to obtain insights into antibiotic resistance genes at the whole genome level.However, few studies have been conducted to comparatively analyze V. parahaemolyticus genomes of aquatic animal origins [16][17][18][19][20].
In our previous research, a number of V. parahaemolyticus strains were collected and analyzed from various aquatic animal species [21].Among these, V. parahaemolyticus B2-28, N9-20, and N2-5 isolates exhibited MDR phenotypes.Hence, we asked what the genome features of the isolates could be and what the molecular mechanism underlying the resistance phenotypes could be.Thus, the primary objectives of this study were (1) to determine draft genome sequences of the three V. parahaemolyticus isolates of aquatic animal origins; (2) to identify MGEs as well as virulence-and resistance-related genes in the V. parahaemolyticus genomes; (3) to investigate the survival of the V. parahaemolyticus strains under different concentrations of antibiotics (AMP, KAN, and STR); and (4) to decipher the molecular mechanism underlying the tolerance of the V. parahaemolyticus strains to the antibiotics through biochemistry, comparative genomics, and transcriptomics analyses (Figure 1).The findings of this study will contribute to the genome data of V. parahaemolyticus and facilitate the control of this foodborne pathogen worldwide.V. parahaemolyticus B2-28, N9-20, and N2-5 strains were isolated from two species of shellfish, namely Ruditapes philippinarum and Keenocardium californiense, respectively, and one species of shrimp, Oratosquilla oratori, respectively [21] (Table S1).V. parahaemolyticus strains were routinely incubated in Tryptic Soy Broth (TSB) medium (3% NaCl, pH 8.5) at 37 °C with shaking at 180 rpm.Bacterial growth was measured using the Bioscreen C au-
2.1.2.Genomic DNA Preparation, Sequencing, Assembly, and Annotation At the middle logarithmic growth phase (mid-LGP), V. parahaemolyticus strains were collected by centrifugation at 8000× g for 1 min.Genomic DNA was extracted using the MiniBEST DNA extraction kit (Japan TaKaRa BIO, Dalian, China).The extracted DNA samples were analyzed using agarose gel electrophoresis, and the DNA concentration and purity were measured [22].Sequencing of only high-quality DNA samples (a 260/280 nm absorbance ratio of 1.8-2.0)was conducted by Shanghai Majorbio Bio-pharm Technology Co., Ltd., Shanghai, China utilizing the Illumina HiSeq×10 sequencing platform (Illumina, Santiago, CA, USA).The PE150 (pair-end) sequencing (insert size: 400 bp) was performed with a read length of 150 bp.Three independently prepared DNA samples were used for each of the V. parahaemolyticus isolates.The positive and negative controls were routinely included in the sequencing run by Shanghai Majorbio Bio-pharm Technology Co., Ltd.(Shanghai, China).
The functional attribution of the predicted coding sequences (CDSs) was inferred using the independent Basic Local Alignment Search Tool (BLAST) (https://www.ncbi.nlm.nih.gov/BLAST, accessed on 16 March 2022) against the NCBI (https://www.ncbi.nlm.nih.gov,accessed on 16 March 2022) database of non-redundant proteins and the database of homologous groups (COG) [26].In cases where the CDS did not exhibit a match with any COG function, it was designated as an unknown protein.The programs were executed using the default parameters.

Serotyping and Multi-Locus Sequence Typing (MLST) Analyses
Serotyping and MLST analyses of the V. parahaemolyticus strains were carried out according to the methods described by Yu et al. [31] and González-Escalona et al. [32], respectively.

Phylogenetic Tree Analysis
The complete genome sequences of 64 V. parahaemolyticus strains were downloaded from the GenBank database.Amino acid data sets of single-copy orthologs present in all the V. parahaemolyticus genomes (n = 67) were inferred and aligned using OrthoFinder (version 2.2.6) [33].Additionally, the FastTree software (version 2.1.11)[34] was used to construct a phylogenetic tree using the same method and parameters described in our previous report [35].

Determination of Minimum Inhibitory Concentrations (MICs) of Antibiotics
V. parahaemolyticus isolates were individually subjected to broth dilution testing (microdilution) following the guidelines by the Clinical and Laboratory Standards Institute (CLSI, M2-A9, 2006, Berwyn, PA, USA).The MICs of antibiotics against the V. parahaemolyticus isolates were determined, including the AMP, KAN, and STR (Sinopharm Chemical Reagent Co., Ltd., Shanghai, China).The MIC is defined as the concentration of a drug that inhibits the visible growth of the bacterium being investigated under defined test conditions.Escherichia coli K-12 (Institute of Industrial Microbiology, Shanghai, China) was used as a quality control strain [21].
Finally, the samples were observed using a thermal field emission SEM (Hitachi, SU5000, Tokyo, Japan) with accelerating voltages ranging from 5 to 10 kV.

Illumina RNA Sequencing
A final concentration of AMP (512 µg/mL), KAN (64 µg/mL), or STR (16 µg/mL) was added into V. parahaemolyticus B2-28, N9-20, and N2-5 culture at the mid-LGP, respectively, and then incubated at 37 • C for 2 h.The bacterial cells were then collected by centrifugation at 12,000 rpm for 2 min.Then, the RNA extraction and quality control processes, as well as the construction of sequencing libraries, were conducted by Shanghai Majorbio Biopharm Technology Co., Ltd. in Shanghai, China using the Illumina HiSeq 2500 platform (Illumina, Santiago, CA, USA).To ensure data accuracy, three replicates were performed for every sample.
Gene expression was determined using the RNA-Seq analysis tool called Expectation-Maximization (RSEM, accessible at http://deweylab.github.io/RSEM/,accessed on 11 December 2022).Genes meeting the criteria of fold changes ≥2.0 or ≤0.5 and p-values < 0.05 were classified as differentially expressed genes (DEGs) compared to the untreated control.A gene set enrichment analysis (GSEA) was performed if the enrichment test p-values were less than 0.05 [38].
The expression of representative DEGs was analyzed using real-time reverse transcription PCR (RT-qPCR) assay [38].Oligonucleotide primers for the RT-qPCR assay were synthesized by Sangon Biotech (Shanghai) Co., Ltd. in Shanghai, China (Table S7).

Data Analysis
The SPSS software (version 22, IBM, Armonk, NY, USA) was used to analyze the data.A one-way analysis of variance was utilized to compare means and sample changes by employing the least-significant difference (LSD) method and homogeneity of variance test with a significance level of p < 0.05.All experiments in this study were performed in triplicate.
The STs of V. parahaemolyticus N9-20, N2-5, and B2-28 isolates were determined by an MLST analysis, and the results show that the V. parahaemolyticus N9-20 and N2-5 isolates belong to ST-1817 and ST-2192, respectively.Notably, the V. parahaemolyticus B2-28 isolate was a new ST type that has not ever been reported.

The Genome Features of the Three V. parahaemolyticus Isolates of Aquatic Animal Origins 3.2.1. General Genome Features
The average nucleotide identity (ANI) values of the three V. parahaemolyticus genomes varied from 98.43% to 98.63%, all surpassing the threshold (94-96%) required for species identification [41].The draft genome sequences of the three V. parahaemolyticus isolates were determined using the Illumina Hiseq×10 sequencing platform (Figure 3).Approximately 173,337-494,024 clean single reads were obtained.The final assembly yielded 42-324 scaffolds with an average sequencing depth ranging from 235.08-fold to 302.48-fold.The obtained genome sizes ranged from 5.0 Mb to 5.4 Mb with average GC contents ranging from 45.14 to 45.40%.A total of 4709-5610 protein-coding genes were predicted.Among these, approximately 3886-4626 genes were classified into 21 functional catalogues against the COG database.Remarkably, the genes with unknown functions accounted for the highest percentages (17.48-21.17%)(Table 1, Figure 3).The obtained genome sizes ranged from 5.0 Mb to 5.4 Mb with average GC contents ranging from 45.14 to 45.40%.A total of 4709-5610 protein-coding genes were predicted.Among these, approximately 3886-4626 genes were classified into 21 functional catalogues against the COG database.Remarkably, the genes with unknown functions accounted for the highest percentages (17.48-21.17%)(Table 1, Figure 3).
The unique 17-mers within the sequencing data of the V. parahaemolyticus N9-20 and N2-5 isolates displayed a clear single peak in frequency.This peak followed a typical Poisson distribution, indicating a lower presence of repetitive DNA in the genomes of V. parahaemolyticus N9-20 and N2-5.Conversely, the V. parahaemolyticus B2-28 genome had a certain degree of heterozygosity with the observed taper at the end of the peak, which suggested relatively more repetitive DNA in its genomes (n = 183) (Figure S1, Table S2).Multiple repetitions were detected towards the terminus of scaffolds (n = 20-183, <1.5 Kb) (Table S2), suggesting that the genome assembly was not fully comprehensive and included numerous gaps.As a result of the limitations of second-generation Illumina short-read sequencing, the unaligned regions between scaffolds consist of recurring sequences [42].The three V. parahaemolyticus genomes contained various MGEs, including GIs (n = 4-14), prophage gene clusters (n = 0-1), Ins (n = 4-27), and ISs (n = 1-2).This suggests the possibility of HGT, which mediates DNA transmission among the bacterial population, facilitated by these MGEs during the evolution of V. parahaemolyticus genomes.It is worth noting that only one of the identified MGEs (IS002 in V. parahaemolyticus B2-28 genome) was found at the end of one of the scaffolds (Tables S3 and S4), indicating that the assembly of draft genomes did not result in the absence of the identified MGEs (except for IS002).

Putative MGEs Genomic Islands
GIs are instrumental in shaping the evolutionary trajectory of V. parahaemolyticus genomes as they acquire new biological characteristics via HGT [43].In this study, 24 GIs (3808-46,379 bp) were identified in the three V. parahaemolyticus genomes, each of which contained 4-14 GIs, carrying 7-49 genes (Figure S2, Table S3).
Interestingly, the genes encoding various functional proteins were found within GIs in three V. parahaemolyticus genomes, e.g., hydrolysis enzymes, chemotaxis proteins, cold shock proteins, stress regulators, transposases, and resistance-related proteins.For example, in the genome of V. parahaemolyticus N2-5, a cold shock protein (Vp_N2-5_0954) was identified in GI3.This protein acts as a global regulator of gene expression and is involved in bacterial growth under various conditions, including adaptation to stress and response to virulence [44].
In 2002, Chamblee et al. sequenced and annotated the complete phage Mu genome that is a temperate phage of a rather wide host range among enteric bacteria [49].Recently, Xu et al. reported three prophages, namely Vibrio phage K139, Pseudomonas phage D3, and Vibrio phage fs2, in the V. parahaemolyticus L7-7, N1-22, N3-33, N4-46, N8-42, and Q8-15 strains isolated from six species of edible aquatic animals [35].In this study, our results coupled with those of a previous study [35] that demonstrated prophage diversity in V. parahaemolyticus genomes.V. parahaemolyticus may have acquired phages from different genera via HGT.Notably, in the genome of V. parahaemolyticus N2-5, ABC transporters were derived from the Enterobacteria phage Mu homologue.These transporter proteins utilize ATP to facilitate the absorption of a wide range of molecules, including nutrients, drugs, antibiotics, and a variety of other compounds [50].

Integrons
Ins enable bacteria to acquire, store, and interchange antibiotic resistance genes [51].Mobile Ins were widely observed in environments where there was a protracted interaction with selective factors like detergents, antibiotics, and heavy metals [52].In this study, all three V. parahaemolyticus genomes contained Ins (n = 4-27) ranging from 662 bp to 227,599 bp.Among them, there were only two complete Ins and a number of gene cassettes (n = 36) (Figure S4, Table S3).
The V. parahaemolyticus B2-28 genome had two ISs, namely an IS3 (1228 bp) encoding two IS3 family transposases and a partial IS5 (962 bp) encoding an AraC family transcriptional regulator.The latter was also found in the V. parahaemolyticus N9-20 genomes.
Taken together, the identified MGEs (n = 69) in this study carrying many genes may have been important drivers of genome evolution and speciation in V. parahaemolyticus.

Putative Virulence-Associated Genes
Previous studies have indicated that V. parahaemolyticus isolates that do not have the tdh and/or trh genes also exhibit significant cytotoxicity towards human gastrointestinal cells, suggesting the presence of alternative virulence-associated genes [61].In this study, several putative genes related to virulence (n = 43-45) were identified in the three V. parahaemolyticus genomes (Table S6).
The BLAST search for heavy metal resistance-associated genes revealed that V. parahaemolyticus B2-28, N9-20, and N2-5 contained the fecE and pfr genes, which were responsible for the tolerance to nickel (Ni), copper (Cu), manganese (Mn), iron (Fe), and cobalt (Co).Notably, the PA0320 gene (Vp_N2_5_0188), which is responsible for the tolerance to Cd 2+ and Hg 2+ , was identified in the V. parahaemolyticus N2-5 genome (Table 2).The gene (Vp_N9-20_4075, GI 6) encoding a short-chain dehydrogenase/reductase SDR family member was identified in V. parahaemolyticus N9-20, which was involved in the response to Cd 2+ stress in Pleurotus eryngii [74], which is consistent with the resistance phenotype to the Cd 2+ of V. parahaemolyticus N9-20.
Variations in resistance genes, genetic diversity, and environmental pressures [75], can lead to distinctions between resistance phenotypes and genotypes.The MIC values of the antibiotics against the V. parahaemolyticus B2-28, N9-20, and N2-5 isolates were determined (Table 3).Remarkably, the observed MICs of V. parahaemolyticus B2-28 against AMP and STR were 100,000 µg/mL and 128 µg/mL, respectively; the MIC values of V. parahaemolyticus N9-20 against AMP, KAN, and STR were 50,000 µg/mL, 256 µg/mL, and 128 µg/mL, respectively; the MIC values of V. parahaemolyticus N2-5 against AMP, KAN, and STR were 50,000 µg/mL, 128 µg/mL, and 128 µg/mL, respectively.The growth curves of the V. parahaemolyticus B2-28, N9-20, and N2-5 isolates were determined at different concentrations of AMP, KAN, and STR in the TSB medium (3% NaCl, pH 8.5) at 37 • C (Figure 4).The growth of V. parahaemolyticus B2-28 was completely inhibited at a concentration of 100,000 µg/mL of AMP.The growth of bacteria was inhibited at concentrations of 50,000 µg/mL and 2048 µg/mL of AMP, resulting in extended lag phases of 10 h and 4 h, respectively.The bacterial biomass reached its peak with OD 600 values of 1.12 and 1.22 after 44 h and 40 h, respectively.Only a minor decline in growth was observed at 512 µg/mL of AMP when compared to the control group (0 µg/mL of AMP) (Figure 4A).Under the AMP (512 µg/mL) condition, the observed fatality rate of V. parahaemolyticus B2-28 was 10.38% (Table 3).At a concentration of KAN of 256 µg/mL, the complete inhibition of V. parahaemolyticus N9-20 was observed (Figure 4B).The bacteria displayed a delayed growth pattern at 128 µg/mL and 64 µg/mL of KAN, with the lag phase being extended to 14 h and 6 h, respectively.Furthermore, the maximum OD600 values were observed at 48 h, reaching 0.91 and 1.25 at 128 µg/mL and 64 µg/mL of KAN, respectively.Only a slight decrease in growth was observed at 64 µg/mL of KAN compared to the control group (0 µg/mL of KAN).Under the KAN (64 µg/mL) condition, V. parahaemolyticus N9-20 showed a fatality rate of 10.02% (Table 3).
The growth of V. parahaemolyticus N2-5 was completely inhibited at a concentration of 128 µg/mL of STR (Figure 4C).Bacterial growth was inhibited at 64 µg/mL and 32 µg/mL of STR with extended lag phases of 12 h and 6 h, respectively.Additionally, the bacterial biomass reached its maximum with OD600 values of 1.10 and 1.25 at 48 h and 30 h, respectively.Growth only slightly decreased at 16 µg/mL of STR compared to the control (0 µg/mL STR).Under the STR (16 µg/mL) condition, the observed mortality rate of V. parahaemolyticus N2-5 was 10.45% (Table 3).

The Cell Membrane Permeability and Fluid and Cell Surface Hydrology of the V. parahaemolyticus Isolates under the Sub-LCs of Antibiotics
External Cell Membrane Permeability Based on the above results, the effects of the sub-LCs of AMP (512 µg/mL), KAN (64 µg/mL), and STR (16 µg/mL) were further determined on the CMP, CMF, and CSH of the V. parahaemolyticus isolates.
The permeation of antibiotics through the outer membrane of Gram-negative bacteria is facilitated by porin channels [82].Specifically, hydrophilic molecules like β-lactams, TET, and some fluoroquinolones are highly influenced by the changes in the permeability of the outer membrane, as they rely on porins to cross this barrier [83].In this study, Nphenyl-1-naphthylamine was used as a probe to detect the external cell membrane permeability.
As shown in Figure 5A, after being treated with AMP (512 µg/mL) for 1.5 h, a consistent decrease in fluorescence intensity was observed.The ECMP of V. parahaemolyticus B2-28 was significantly increased by 1.04-to 1.09-fold when compared with the control group (p < 0.05).At a concentration of KAN of 256 µg/mL, the complete inhibition of V. parahaemolyticus N9-20 was observed (Figure 4B).The bacteria displayed a delayed growth pattern at 128 µg/mL and 64 µg/mL of KAN, with the lag phase being extended to 14 h and 6 h, respectively.Furthermore, the maximum OD 600 values were observed at 48 h, reaching 0.91 and 1.25 at 128 µg/mL and 64 µg/mL of KAN, respectively.Only a slight decrease in growth was observed at 64 µg/mL of KAN compared to the control group (0 µg/mL of KAN).Under the KAN (64 µg/mL) condition, V. parahaemolyticus N9-20 showed a fatality rate of 10.02% (Table 3).
The growth of V. parahaemolyticus N2-5 was completely inhibited at a concentration of 128 µg/mL of STR (Figure 4C).Bacterial growth was inhibited at 64 µg/mL and 32 µg/mL of STR with extended lag phases of 12 h and 6 h, respectively.Additionally, the bacterial biomass reached its maximum with OD 600 values of 1.10 and 1.25 at 48 h and 30 h, respectively.Growth only slightly decreased at 16 µg/mL of STR compared to the control (0 µg/mL STR).Under the STR (16 µg/mL) condition, the observed mortality rate of V. parahaemolyticus N2-5 was 10.45% (Table 3).

The Tolerance Mechanisms of V. parahaemolyticus Isolates from Aquatic Animals to the
Sub-LCs of the Three Antibiotics 3.4.1.The Cell Membrane Permeability and Fluid and Cell Surface Hydrology of the V. parahaemolyticus Isolates under the Sub-LCs of Antibiotics External Cell Membrane Permeability Based on the above results, the effects of the sub-LCs of AMP (512 µg/mL), KAN (64 µg/mL), and STR (16 µg/mL) were further determined on the CMP, CMF, and CSH of the V. parahaemolyticus isolates.
The permeation of antibiotics through the outer membrane of Gram-negative bacteria is facilitated by porin channels [82].Specifically, hydrophilic molecules like β-lactams, TET, and some fluoroquinolones are highly influenced by the changes in the permeability of the outer membrane, as they rely on porins to cross this barrier [83].In this study, N-phenyl-1naphthylamine was used as a probe to detect the external cell membrane permeability.
As shown in Figure 5A, after being treated with AMP (512 µg/mL) for 1.5 h, a consistent decrease in fluorescence intensity was observed.The ECMP of V. parahaemolyticus B2-28 was significantly increased by 1.04-to 1.09-fold when compared with the control group (p < 0.05).After the treatment of KAN (64 µg/mL), the fluorescence intensity consistently decreased over a period of 4 h.The ECMP of V. parahaemolyticus N9-20 was significantly increased by 1.01-and 1.12-fold at 2 h and 3.5 h compared with the control group (p < 0.05) (Figure 5B).
Similarly, the fluorescence intensity showed an overall downward trend after the treatment with STR (16 µg/mL) for 0 h to 4 h.The ECMP of V. parahaemolyticus N2-5 was significantly increased by 1.05-to 1.16-fold compared with the control group (Figure 5C) (p < 0.05).

Internal Cell Membrane Permeability
In this study, the O-nitrophenyl-β-D galactopyranoside was used as a probe to monitor internal cell membrane permeability.As shown in Figure 5D, no significant difference in the ICMP of V. parahaemolyticus B2-28 was observed after being treated with AMP (512 µg/mL) for 0 h to 2.0 h and for 3.0 h to 4.0 h compared with the control group (p > 0.05).However, its ICMP increased by 18.60-fold at 2.5 h (p < 0.05).
Similarly, no significant difference in the ICMP of V. parahaemolyticus N9-20 was observed after being treated with KAN (64 µg/mL) for 0 h to 1.0 h compared with the control group (p > 0.05) (Figure 5E).However, during the longer treatment time (1.5 to 3.0 h), its CIMP increased by 2.95-to 4.22-fold (p < 0.05).
Likewise, after being treated with STR (16 µg/mL) for 0 h to 1.0 h and for 2.5 h to 4.0 h, the ICMP of V. parahaemolyticus N2-5 showed no significant change compared with the control group (p < 0.05) (Figure 5F).However, at 1.5 to 2.0 h, its ICMP increased by 3.13and 4.87-fold (p < 0.05).

Cell Surface Hydrology and Cell Membrane Fluid
The CSH plays an important role in bacterial adhesion to abiotic and biotic surfaces as well as penetration into host tissues [84].As shown in Figure 6A, when compared with After the treatment of KAN (64 µg/mL), the fluorescence intensity consistently decreased over a period of 4 h.The ECMP of V. parahaemolyticus N9-20 was significantly increased by 1.01-and 1.12-fold at 2 h and 3.5 h compared with the control group (p < 0.05) (Figure 5B).
Similarly, the fluorescence intensity showed an overall downward trend after the treatment with STR (16 µg/mL) for 0 h to 4 h.The ECMP of V. parahaemolyticus N2-5 was significantly increased by 1.05-to 1.16-fold compared with the control group (Figure 5C) (p < 0.05).

Internal Cell Membrane Permeability
In this study, the O-nitrophenyl-β-D galactopyranoside was used as a probe to monitor internal cell membrane permeability.As shown in Figure 5D, no significant difference in the ICMP of V. parahaemolyticus B2-28 was observed after being treated with AMP (512 µg/mL) for 0 h to 2.0 h and for 3.0 h to 4.0 h compared with the control group (p > 0.05).However, its ICMP increased by 18.60-fold at 2.5 h (p < 0.05).
Similarly, no significant difference in the ICMP of V. parahaemolyticus N9-20 was observed after being treated with KAN (64 µg/mL) for 0 h to 1.0 h compared with the control group (p > 0.05) (Figure 5E).However, during the longer treatment time (1.5 to 3.0 h), its CIMP increased by 2.95-to 4.22-fold (p < 0.05).
Likewise, after being treated with STR (16 µg/mL) for 0 h to 1.0 h and for 2.5 h to 4.0 h, the ICMP of V. parahaemolyticus N2-5 showed no significant change compared with the control group (p < 0.05) (Figure 5F).However, at 1.5 to 2.0 h, its ICMP increased by 3.13and 4.87-fold (p < 0.05).

Cell Surface Hydrology and Cell Membrane Fluid
The CSH plays an important role in bacterial adhesion to abiotic and biotic surfaces as well as penetration into host tissues [84].As shown in Figure 6A, when compared with the control group, the CSH of V. parahaemolyticus B2-28 was significantly increased by 4.20-fold after the treatment with AMP (512 µg/mL) for 2 h (p < 0.001).In contrast, the treatment with KAN (64 µg/mL) for 2 h significantly decreased the CSH of V. parahaemolyticus N9-20 by 1.89-fold compared with the control group (p < 0.001).Similarly, the treatment with STR (16 µg/mL) for 2 h significantly reduced the CSH of V. parahaemolyticus N2-5 by 1.53-fold compared to the control group (p < 0.05).The CMF plays a key role in the action of membrane-active antibiotics [85].As shown in Figure 6B, when compared with the control group, the CMF of V. parahaemolyticus N9-20 was significantly increased by 1.12-fold (p < 0.001) after the treatment with KAN (64 µg/mL) for 2 h.The treatment with STR (16 µg/mL) for 2 h significantly increased the CMF of V. parahaemolyticus N2-5 by 1.19-fold (p < 0.05).In contrast, the treatment with AMP (512 µg/mL) for 2 h significantly decreased the CMF of V. parahaemolyticus B2-28 by 1.07-fold (p < 0.05).
Taken together, under the sub-LCs of antibiotics, the ECMP, ICMP, CSH, and CMF of the V. parahaemolyticus B2-28, N9-20, and N2-5 isolates were significantly altered (p < 0.05), suggesting that there was a possible change in the bacterial cellular structure.

The Cell Structure Change of the Three V. parahaemolyticus Isolates under the Sub-LCs of Antibiotics
As shown in Figure 7, the control groups displayed intact rod cells with a flat and transparent surface structure, whereas the treatment groups exhibited a slight contraction.Some of the cells in the treatment groups appeared to be deformed, showing folds on their surfaces.For example, when compared with the control groups, after the treatment with AMP (512 µg/mL) for 2 h, the cell surface of V. parahaemolyticus B2-28 shrunk slightly and was subtly depressed (Figure 7A); for V. parahaemolyticus N9-20, the treatment for 2 h using KAN (64 µg/mL) resulted in the slight shrinking in the bacterial cell surface (Figure 7B).A similar case was observed for V. parahaemolyticus N2-5 after the treatment with STR (16 µg/mL) for 2 h (Figure 7C).The CMF plays a key role in the action of membrane-active antibiotics [85].As shown in Figure 6B, when compared with the control group, the CMF of V. parahaemolyticus N9-20 was significantly increased by 1.12-fold (p < 0.001) after the treatment with KAN (64 µg/mL) for 2 h.The treatment with STR (16 µg/mL) for 2 h significantly increased the CMF of V. parahaemolyticus N2-5 by 1.19-fold (p < 0.05).In contrast, the treatment with AMP (512 µg/mL) for 2 h significantly decreased the CMF of V. parahaemolyticus B2-28 by 1.07-fold (p < 0.05).
Taken together, under the sub-LCs of antibiotics, the ECMP, ICMP, CSH, and CMF of the V. parahaemolyticus B2-28, N9-20, and N2-5 isolates were significantly altered (p < 0.05), suggesting that there was a possible change in the bacterial cellular structure.

The Cell Structure Change of the Three V. parahaemolyticus Isolates under the Sub-LCs of Antibiotics
As shown in Figure 7, the control groups displayed intact rod cells with a flat and transparent surface structure, whereas the treatment groups exhibited a slight contraction.Some of the cells in the treatment groups appeared to be deformed, showing folds on their surfaces.For example, when compared with the control groups, after the treatment with AMP (512 µg/mL) for 2 h, the cell surface of V. parahaemolyticus B2-28 shrunk slightly and was subtly depressed (Figure 7A); for V. parahaemolyticus N9-20, the treatment for 2 h using KAN (64 µg/mL) resulted in the slight shrinking in the bacterial cell surface (Figure 7B).A similar case was observed for V. parahaemolyticus N2-5 after the treatment with STR (16 µg/mL) for 2 h (Figure 7C).

The Differential Transcriptomes of the Three V. parahaemolyticus Isolates Induced by the Sub-LCs of Antibiotics
Based on the obtained results, V. parahaemolyticus B2-28, N9-20, and N2-5 grown at mid-LGP in the TSB medium at 37 °C were treated with AMP (512 µg/mL), KAN (64 µg/mL), and STR (16 µg/mL) for 2 h, respectively, and gene expression changes at the global genome level were investigated using Illumina HiSeq 2500 sequencing technology.
The DEGs (n = 11) that were related to the FAS II pathway [87,88] in the fatty acid biosynthesis were all significantly up-regulated (2.018-to 2.803-fold) (p < 0.05).Of these, the FabF protein (Vp_B2_28_0092) was extensively proven as the desired focus in Grampositive bacteria targeted by natural compounds such as platensimycin, platencin, and fasamycin A and B. [87].The accC (Vp_B2_28_4804) and accA (Vp_B2_28_1017) genes encoded an ACCase transferase subunit alpha and an ACCase biotin carboxylase subunit, respectively.The ACCase subunit is highly conserved in Gram-positive and Gram-negative bacteria, suggesting that it is possible to find ACCase inhibitors with broad-spectrum antibacterial activity [88].Moreover, the necessity of ACCase activity for bacterial growth and survival has been clearly established [88].
AMP stress also triggered significant changes in the other three metabolic pathways in V. parahaemolyticus B2-28.All DEGs (n = 14) were significantly up-regulated (2.028to 12.837-fold) in the GSTM (p < 0.05).For instance, the ectABC genes (Vp_B2_28_3088, Vp_B2_28_3089, and Vp_B2_28_3090) encoding a diaminobutyrate acetyltransferase were diaminobutyrate-2-oxoglutarate transaminase genes, and an ectoine synthase in the ectoine biosynthesis in the GSTM was significantly up-regulated (2.028-to 2.065-fold) (p < 0.05).Ectoines, which are extensively employed by microorganisms, are noteworthy compatible solutes.They effectively safeguard protein functionality when faced with diverse chal-lenges, impact membrane fluidity, and stabilize lipid bilayers [90].This was consistent with the increased CMF of the V. parahaemolyticus isolates treated by AMP in this study.
Taken together, under the sub-LC of AMP stress, V. parahaemolyticus B2-28 significantly up-regulated the enzymatic activity of PBPs to eliminate AMP's binding ability and enhanced the assimilation of sulfate, fructose transport, ectoines, and endogenous cadaverine biosynthesis to reduce AMP damage to cells.
In contrast, most DEGs (n = 22 of 25) were significantly up-regulated in amino acid metabolism (2.010 to 7.057-fold) (p < 0.05), such as the hisABCDFGH(IE) (Vp_N9_20_1730 to Vp_N9_20_1737) and leuABCD (Vp_N9_20_2872 to Vp_N9_20_2875) genes, which were related to L-histidine and L-leucine synthesis, respectively.Also, the cysEK genes (Vp_N9_20_4215 and Vp_N9_20_1922) [96,97] encoding a serine O-acetyltransferase and a PLP-dependent cysteine synthase family protein related to cysteine were significantly up-regulated (2.222 to 2.508-fold) (p < 0.05).The preservation of bacterial cells from redox harm and the development of antibiotic resistance (such as amphotericin B, KAN, and chloramphenicol) are closely tied to cysteine and its associated metabolites [97].
In biofilm formation, six DEGs encoding the type VI secretion system (T6SS) were significantly down-regulated (0.289 to 0.5-fold) (p < 0.05).The T6SS of V. cholerae and Pseudomonas aeruginosa has been shown to be involved in the export of hemolysin-coregulated proteins and valine-glycine repeat proteins [54].In this study, the T3SS transcriptional regulator ExsA (Vp_N9_20_2745) was significantly up-regulated by 7.082-fold (p < 0.05).The majority of T3SS1 genes require the master regulator ExsA for their expression [98].
Taken together, V. parahaemolyticus N9-20 employed multiple strategies to cope with the sub-LC of KAN (64 µg/mL) stress: (1) it reduced soluble ligand PPBPs and the transmembrane transport of choline, betaine, tungstate, L-arabinose, and ribose; and (2) it decreased carbon source utilization to increase the CMP; (3) it promoted the accumulation of amino acid metabolites, such as L-histidine, L-leucine, glutathione, and cystine, to improve cell survivability.
The Major Changed Metabolic Pathways in V. parahaemolyticus N2-5 under STR Stress Approximately 5.6% (293 of 5191 genes) of the bacterial genes were differentially expressed in V. parahaemolyticus N2-5 after being treated by the sub-LC of STR (16 µg/mL) for 2 h compared to the control group.Of these, 90 DEGs showed lower transcriptional levels (fold change ≤ 0.5), while 203 DEGs were up-regulated (fold change ≥ 2.0) (Figure 8C,F).Approximately six significantly altered metabolic pathways were identified in V. parahaemolyticus N2-5, including ABC transporters, GSTM, lysine biosynthesis, the MAPK signaling pathway, sulfur metabolism, and nitrogen metabolism (Table S10).
The DEGs (n = 3) in the lysine biosynthesis were also up-regulated (2.156-to 3.935-fold) (p < 0.05), such as the lysC gene (Vp_N2_5_4367) encoding a Lysine-sensitive aspartokinase 3, and the thrA gene (Vp_N2_5_3440) encoding a bifunctional aspartate kinase/homoserine dehydrogenase I. Lysine plays a crucial role in the adaptation and tolerance to environmental stresses in diverse organisms [99].
In terms of energy metabolism, the DEGs (n = 7) involved in sulfur metabolism and nitrogen metabolism were significantly up-regulated (2.03-to 3.044-fold) (p < 0.05), such as the metA gene (Vp_N2_5_3773) encoding a homoserine O-succinyltransferase, the fccA gene (Vp_N2_5_2651) encoding a C-type cytochrome, and the gltBD genes (Vp_N2_5_3431 and Vp_N2_5_3428) encoding a glutamate synthase large and small subunits.
Taken together, V. parahaemolyticus N2-5 employed multiple strategies to cope with the sub-LC of STR (16 µg/mL) stress: (1) it reduced the membrane transport of arginine and the biosynthesis of ectoine (2) and promoted sugar, choline, ribose, and phosphate membrane transport and energy metabolism to reduce STR damage to cells.
Additionally, the expressions of representative DEGs were examined using an RT-PCR assay, and the resulting data were consistent with the transcriptomic analysis (Table S11).
3.4.4.The Molecular Basis Underlying the Antibiotic Tolerance of the V. parahaemolyticus Isolates at the Sub-LCs of the Antibiotics A comparative transcriptomic analysis revealed nine significantly altered metabolic pathways in V. parahaemolyticus B2-28 under the sub-LC of AMP (512 µg/mL).Remarkably, the DEG (blaCARB-17) encoding carbenicillin-hydrolyzing class A beta-lactamase was highly up-regulated by 6.409-fold (p < 0.05).Interestingly, PBPs (mrcAB and mrdA) in peptidoglycan biosynthesis were significantly up-regulated by 2.003-to 2.472-fold (p < 0.05).In Gram-negative bacteria, the most prevalent mechanism for the development of β-lactam resistance is the production of β-lactamases, followed by altered permeability, the extrusion of efflux pumps, and altered PBPs [100], which is consistent with the mechanism of the β-lactam resistance of V. parahaemolyticus B2-28 in this study.In sulfur metabolism, DEGs participate in the assimilation of sulfate (cysCDNHIJ), and tetrathionate reductase (ttrBCA) was significantly up-regulated (p < 0.05).The DEGs in fructose and mannose metabolism and PTS (fruABK), the ectoine biosynthesis (ectABC) in the GSTM, and the endogenous cadaverine (ldcC, cadA) in the TPPAB were also significantly up-regulated (p < 0.05).These results indicated the increased substance synthesis for energy conservation and stringent response regulation in V. parahaemolyticus B2-28 under AMP (512 µg/mL) stress.
Interestingly, DEGs (fabAFGHV, fadD, and accACD) that participate in the FAS II pathway in the fatty acid biosynthesis were up-regulated in V. parahaemolyticus B2-28 under AMP stress (p < 0.05).These results suggest that AMP may target the FAS II pathway and thus affect cell membrane synthesis.
A comparative transcriptomic analysis revealed nine significantly altered metabolic pathways in V. parahaemolyticus N9-20 under the sub-LC of KAN (64 µg/mL).Remarkably, the DEGs encoding the molecular chaperone DnaK (dnaK), DnaJ (dnaJ), HtpG (htpG), and the Hsp20 family protein (ibpA) were greatly up-regulated by 15.337-, 16.871-, 16.065, and 87.874-fold (p < 0.05), respectively.KAN, a type of aminoglycoside antibiotic, disrupts the process of protein synthesis by binding to the 16S rRNA located in the A site of the 30S ribosomal subunit.This interaction hinders the accurate selection of cognate tRNAs, leading to the synthesis of abnormal proteins [101].Zhang et al. reported that stress-related proteins (GroS, DnaK, GroL, HtpG, ClpB, HslU, and DnaJ) are hub proteins that significantly increase to reduce the pressure from the misreading of mRNA caused by KAN [102].In this study, interestingly, hisABCDFGH(IE) in the histidine metabolism, leuABCD in the VLIB, and gshAB and cysEK in the CMM were significantly up-regulated (p < 0.05).The sucCD in the C5-BDAM and PHA (phbBC and atoB) in the butanoate metabolism were significantly regulated (p < 0.05).These results indicate that V. parahaemolyticus N9-20 utilized PHA for carbon storage and promoted the accumulation of amino acid metabolites under KAN (64 µg/mL) stress.
The DEGs encoding a homoserine O-succinyltransferase (metA) and a c-type cytochrome (fccA) in the sulfur metabolism and the large (gltB) and small (gltD) subunits of glutamate synthase in the nitrogen metabolism were significantly up-regulated by 2.030to 3.044-fold (p < 0.05) in V. parahaemolyticus N2-5.In the ABC transporter, the DEGs encoding arginine (artIMP) were significantly down-regulated (p < 0.05).On the other hand, the DEGs encoding ribose (rbsAD), sugar (malEG, araH), phosphate (pstS and afuB), and choline (proW) were significantly up-regulated (p < 0.05).These results suggest the possible altered energy metabolism for the transport of ribose, sugar, choline, and phosphate in V. parahaemolyticus N2-5 under STR stress.

Foods 2024 , 27 Figure 1 .
Figure 1.A flow chart of the investigation strategy in this study.

2. Materials and Methods 2 . 1 .
The Characterization of Genome Features of the V. parahaemolyticus Isolates 2.1.1.V. parahaemolyticus Strains and Culture Conditions

Figure 1 .
Figure 1.A flow chart of the investigation strategy in this study.

Figure 2 .
Figure 2. Genome-wide homologous single-copy gene evolutionary tree with 67 V. parahaemolyticus strains.From the inner circle to the outer circle: isolation location, isolation time, host information, ST type, and serotype of 67 V. parahaemolyticus strains, respectively.

Figure 2 .
Figure 2. Genome-wide homologous single-copy gene evolutionary tree with 67 V. parahaemolyticus strains.From the inner circle to the outer circle: isolation location, isolation time, host information, ST type, and serotype of 67 V. parahaemolyticus strains, respectively.

Foods 2024 , 27 Figure 3 .
Figure 3. Genome circle maps of the three V. parahaemolyticus isolates.(A,B) represent the larger and smaller chromosomes of the three V. parahaemolyticus isolates, respectively.Circles from the inside to the outside: the GC contents, the GC skew values (positive in purple and negative in green), the reference genome of V. parahaemolyticus RIMD 2210633 (GenBank accession numbers NC_004603.1 and NC_004605.1),the genomes of V. parahaemolyticus B2-28, N9-20, and N2-5, as well as the CDSs (positioned on the positive and negative strands accordingly), respectively.

Figure 3 .
Figure 3. Genome circle maps of the three V. parahaemolyticus isolates.(A,B) represent the larger and smaller chromosomes of the three V. parahaemolyticus isolates, respectively.Circles from the inside to the outside: the GC contents, the GC skew values (positive in purple and negative in green), the reference genome of V. parahaemolyticus RIMD 2210633 (GenBank accession numbers NC_004603.1 and NC_004605.1),the genomes of V. parahaemolyticus B2-28, N9-20, and N2-5, as well as the CDSs (positioned on the positive and negative strands accordingly), respectively.

Foods 2024 ,
13,  x FOR PEER REVIEW 15 of 27 the control group, the CSH of V. parahaemolyticus B2-28 was significantly increased by 4.20fold after the treatment with AMP (512 µg/mL) for 2 h (p < 0.001).In contrast, the treatment with KAN (64 µg/mL) for 2 h significantly decreased the CSH of V. parahaemolyticus N9-20 by 1.89-fold compared with the control group (p < 0.001).Similarly, the treatment with STR (16 µg/mL) for 2 h significantly reduced the CSH of V. parahaemolyticus N2-5 by 1.53fold compared to the control group (p < 0.05).

3. 4 . 3 .
The Differential Transcriptomes of the Three V. parahaemolyticus Isolates Induced by the Sub-LCs of Antibiotics

Table 1 .
General features of the three V. parahaemolyticus genomes.

Table 1 .
General features of the three V. parahaemolyticus genomes.

Table 2 .
The antimicrobial resistance-and heavy metal resistance-associated genes identified in the three V. parahaemolyticus genomes.
3.3.The MICs and Sublethal Concentrations (Sub-LCs) of the Three V. parahaemolyticus Isolates against Antibiotics

Table 3 .
The MIC values and sub-LCs of the V. parahaemolyticus isolates against antibiotics. V.