Identification of the AHPND positive strains
Molecular identification and characterization of suspected AHPND positive V. parahaemolyticus isolates were done using 16S rRNA, ldh, AP3 and AP4 primers PCR (Figure 1). MSR16 (isolated in 2016) and MSR17 (isolated in 2017) strains were finally sequenced for whole genome sequencing.
Features of the assembled genomes
The genomes were assembled into 35 contigs in MSR16 strain and 53 contigs in MSR17 strain. The largest contigs size for MSR16 strain was 1742619 bp; and 1892806 bp for MSR17 strain. The total GC content was 45.19% and 45.09% for MSR16 and MSR17 strains, respectively. The total genome size of MSR16 was 5232129 bp; and 5377931 bp for MSR17. MSR16 was found comprised of two circular chromosome with length of 3,407,823 bp, 1,893,040 bp while genome of MSR17 was comprised with similar two circular chromosome with length of 3,428,841 bp, 1,743,028 bp. Both MSR16 and MSR17 contain a plasmid with a length of 68,345 bp and 66,825 bp, respectively (Figure 2). Comparing the genomes, it was observed that chromosome 2 of MSR16 strain has an extra 100kb region. More information about MSR16 and MSR17 genomes are given in Table 1.
Table 1: Summary of the assembled genomes of two strains (MSR16 & MSR17) of AHPND positive V. parahaemolyticus.
Features
|
VPAHPND MSR16
|
VPAHPND MSR17
|
Contigs
|
108
|
66
|
Largest contigs
|
1742619
|
1892806
|
Total length
|
5232129
|
5377931
|
GC (%)
|
45.19
|
45.09
|
CDS
|
4909
|
4689
|
Gene
|
5090
|
4854
|
tRNA
|
119
|
109
|
misc_RNA
|
51
|
45
|
rRNA
|
10
|
10
|
tmRNA
|
1
|
1
|
The plasmid of MSR16 contains total 87 genes of which 58 genes are hypothetical protein (67%), 5 repeat regions (6%), 7 conjugative transfer proteins (8%), 3 mobile element protein (3%), 2 antirestriction protein (2%), 2 toxin genes (pirA and pirB) and 10 other genes (11%). The plasmid of MSR17 contains total 88 genes of which 57 genes are hypothetical protein (65%), 6 repeat regions (7%), 7 conjugative transfer proteins (8%), 3 mobile element protein (3%), 2 antirestriction protein (2%), 2 toxin genes (pirA and pirB) and 11 other genes (13%).
Out of the RAST server predicted 406 subsystems, MSR16 strain possesses 74 responsible for virulence, disease and defense; five for phages, prophages, transposable elements and plasmids; 28 for iron acquisition and metabolism; and 125 for motility and chemotaxis. While out of the predicted 403 subsystems, MSR17 strain contained 74 responsible for virulence, disease and defense; 10 for phages, prophages, transposable elements and plasmids; 28 for iron acquisition and metabolism; and 119 for motility and chemotaxis (Figure 3). These particular subsystems are the hallmarks for the pathogenicity and both strains were found to have almost similar amount of factors across their genomes. Number of genes associated with the general COG functional categories for both strains are provided in (Figure 4). Both strains are found to possess equivalent number of genes associated with those categories.
MSR16 and MSR17 strains have average nucleotide identity values of 98.57% with V. parahaemolyticus strain M1-1 and 98.65% with V. parahaemolyticus strain 13-306D/4 respectively; they also have an average of 95% ANI values with other AHPND positive strains (Additional file 1). Strains MSR16 and MSR17 have 1403 and 1228 hypothetical genes respectively, whose functional prediction can provide more insights into its pathogenicity and other functional pathways. 144 and 94 unique genes were found in strain MSR16 and MSR17 respectively which are uniquely predicted only for one strain (Additional file 2). MSR17 strain contains unique genes for zona occludens toxin, several transposition proteins, integrase, recombinases etc.; whereas MSR16 strain have genes for several conjugative transfer related proteins, bacteriocin immunity proteins etc. Both strains are predicted to have some exclusive genes for diverse metabolic pathways.
Virulence and antimicrobial resistance genes
Most common 9 virulence factor classes involved in- adherence, antiphagocytosis, enzyme, chemotaxis and motility, iron uptake, quorum sensing, secretion system, toxin, immune envasion were found in the MSR16, while MSR17 possess 8 of these such factors except the factors involved in immune envasion; also few genes in these classes of factors were found absent in these strains (Additional file 3). The thermostable direct haemolysin (tdh), the TDH-related haemolysin (trh) and the two type III secretion systems (T3SS1 and T3SS2) are recognized as major virulence factors in V. parahaemolyticus (17). tdh and trh both genes were not found in these strains but the thermo labile hemolysin (tlh) gene was found. Between two types of T3SS only T3SS1 type was found in MSR16 and MSR17 strain. Both strains possess the plasmid borne pirA and pirB toxins.
Antibiotic resistance genes were predicted against β-lactam, fluoroquinolone, tetracycline, macrolide and cephalosporin antibiotics in MSR16; and MSR17 strain has similar resistance genes except for cephalosporin (Additional file 4). Six and two probable prophage regions were found in MSR16 and MSR17 strains, respectively.
Strains MSR16 and MSR17 have approximately 39 and 27 genomic island (GI) regions respectively (Additional file 5). In strain MSR16, toxin-antitoxin systems like YoeB-YefM, Doc-Phd; antibiotic resistance proteins like FosA (Fosfomycin resistance protein); components of type-I, type-VI secretion systems etc. are found in those genomic islands. Genomic islands of strain MSR17 contains toxin-antitoxin systems like HipA-HipB, YoeB-YefM; type-I, type-III secretion systems; Multidrug resistance efflux pump; several phage and transposon related proteins etc. (Additional file 6).
PathogenFinder tool (18) predicted an overall probability of 0.868 for MSR16 and 0.871 for MSR17 for becoming potential human pathogen, so there is a very high risk of spreading these strains into the human food chain and causing human diseases.
Phylogenetic relationship based on 16S rRNA genes of different AHPND positive V. parahaemolyticus strain
A total of 30 strains were selected for establishing phylogenetic relationship based on 16S rRNA gene sequence (Figure 5). The tree includes 25 V. parahaemolyticus (including MSR16 and MSR17), two V. campbellii and two V. owensii strains which were responsible for AHPND outbreak in the recent years in different regions of the world. V. cholerae was used for outgroup comparison.
The tree illustrates that Vibrio campbellii AB497067 India, Vibrio campbellii AY035896 China strains were located at same cluster. The strains including M13-18718-18H KY229843.1 Australia, NSTH24 KF886635.1 Thailand, NSTH22 KF886633.1 Thailand, NSTH21 KF886632.1 Thailand, L44 KC884620.1 China, CZN-33 KR347269.1 China, MM21 FJ547093.1 China, Vp-D4 MH298544.1 China, and Vp-X 11 MH298549.1 China located at same cluster.
Both of our studied strains (MSR16 and MSR17) located at same cluster and were closely related with one of Indian strain AP1511 indicating that the mutation and evolutionary pattern of MSR16 and MSR17 strains might be analogous to the Indian strain. All the Indian V. parahaemolyticus strains GM171 MG593212.1 India, AP183 MG564744.1 India, SPJ3 MG564747.1 India, G476 MG593229.1 India, GA5114 MG970587.1 India, S24P132 MG762012.1 India, and AP1511 MG564729.1 India located at same cluster. The two Spanish V. parahaemolyticus strains (CECT5305 FM204865 Spain and CECT611 FM204867.1 Spain) separately make a cluster. One of the Chinese V. parahaemolyticus strain 461 JN188418.1 China belonged to an independent lineage.
The strains including Vp-4 MK377081.1 China, Ramsar KJ704113.2 Iran belong to separate cluster. Besides, two AHPND positive V. owensii strains S12A KR086360 India, CHN-25 KR347195 China were located at separate cluster. V. cholerae (msr6) strain was distantly related with our studied strains.
Phylogenetic relationship based on housekeeping genes of different AHPND positive V. parahaemolyticus
A total of 25 strains were selected for establishing phylogenetic relationship based on common housekeeping genes (Figure 6) including (dnaE, dtdS, gyrB, pntA, pyrC, recA, tnaA). The 16S rRNA gene was not included because separate phylogenetic relationship was established based on it. FIM-S1708+ Mexico, FIM-S1392- Mexico and MVP5 Malaysia strains located at same cluster. TUMSAT-DE1-S1 Thailand, NCKU-TN-S02 Thailand, NCKU-TV-5HP Thailand and KS17.S5-1 Malaysia strains located at same cluster. PSU5570 Thailand and ND11 Malaysia strain belong to an independent lineage respectively. The strains M0605 Mexico and TUMSAT-H10-S6 Thailand located at same cluster. NCKU-TV-3HP Thailand, MSR17 Bangladesh, M1-1 Vietnam located at same cluster. MVP3 Malaysia and VP14 India strain located at same cluster. The strains 12-009A/1335 Vietnam, MSR16 Bangladesh located at same cluster. 13-028-A2 Vietnam and NA9 Malaysia strains separately make a cluster. NCKU_CV_CHN_China belong to an independent lineage. TUMSAT-D06-S3 Thailand, 13-306-D4 Mexico and NB09 Malaysia strain located at same cluster.
The phylogenetic tree showed that MSR16 strain was closely related to 12-009A/1335 Vietnam strain which maintain an antibacterial type vi secretion system with versatile effector repertoires (19) suggesting that MSR16 strain is originated from Vietnam. MSR17 strain was closely related to M1-1 Vietnam strain signifying that MSR17 strain might evolves from M1-1 Vietnamese strain. Kumar et al. (2018) reported that M1-1 strain causes a mild form of shrimp acute hepatopancreatic necrosis disease (AHPND) (20). Compared to other virulent strains, the M1-1 genome appeared to express several additional genes, while some genes were missing. These instabilities may be related to the reduced virulence of M1-1.
The tree also illuminates that MSR16 strain arise earlier than MSR17 strain. NA7 Malaysia strain belonged to an independent lineage and distantly related to our studied strains (MSR16 and MSR17) signifying that this strain might be disparate to MSR16 and MSR17.
SNP tree of different AHPND positive V. parahaemolyticus
A total of 37 genome of AHPND positive V. parahaemolyticus strain including MSR16 and MSR17 were selected for establishing SNP relationship (Figure 7). MVP3 Malaysia and MVP5 Malaysia strain located at same cluster. The strains MVP10 Malaysia and 13-028-A2 Vietnam located at same cluster. TUMSAT-H03-S5 Thailand, FIM-S1392 Mexico and TUMSAT-H10-S6 Thailand strain belongs to an independent lineage respectively. The strains NB08 Malaysia and NB10 Malaysia separately make a cluster. The strain including NA8 Malaysia, NA6 Malaysia and NB07 Malaysia located at same cluster. The strain including NB12 Malaysia and NB09 Malaysia located at same cluster.
The strains MSR16 Bangladesh and NA9 Malaysia were closely related and located at same cluster indicating that the mutation and evolutionary pattern of MSR16 might be comparable to Malaysian strains. NA9 strain was extracted from Malaysian aquaculture pond water which causes AHPND in shrimp and impacting Malaysian shrimp aquaculture. While strains M0605 Mexico and MSR17 Bangladesh were closely related and located at same cluster indicating that the mutation and evolutionary pattern of MSR17 might be analogous to the Mexican strain. Gomez-Gil et al. (2014) reported that, several pathogenicity mechanisms were identified on both chromosomes: five iron acquisition systems (hemin, enterobactin, vibrioferrin, and two TonB) and seven secretion systems (two T2SS, one T3SS, two T2/4SS, and two T6SS). At least 14 different toxin genes were annotated, two of which are large repeats in toxin (RTX), as well as nine hemolysins. Gomez-Gil et al. (2014) also detected four plasmid from M0605 strains genome (11).
Strains NA7 Malaysia and NA4 Malaysia located at same cluster. Strain TUMSAT-H03-S5 Thailand strain belonged to an independent lineage. This strains mutation and evolutionary pattern might disparate to MSR16 and MSR17.
ANI (Average Nucleotide Identity) tree of different AHPND positive V. parahaemolyticus strain
A total of 52 genome of AHPND positive V. parahaemolyticus strain including MSR16 and MSR17 were selected for calculating the average nucleotide identity (ANI) (Figure 8). The 11 strains including VP14 India, FIMS1392-, 13-028/A2 Vietnam, TUMSAT-S5 Thailand, TUMAT-S4 Thailand, M0605-Mexico, 13-306-D4 Mexico, MSR17, FIMS1708+ Mexico, TUMSATDE1 Thailand, NCKU-CHN China were located at cluster 2. KS17.S5-1 Malaysia, TUMSATD06 Thailand, 12-297/B Vietnam, 12-009A Vietnam, PSU5579 Thailand and NCKU-S02 China located at cluster same cluster.
The ANI tree clearly illustrates that MSR16 strain belonged to an independent lineage and indicating this strain evolve earlier than MSR17. The reasons for belongs to an independent linage might be the presence of extra 200 kb nucleotide in the genome. The strain MSR17 was closely related to 13-306-D4 Mexico strain signifying that the average nucleotide identity (ANI) of MSR17 is comparable to this Mexican strain.
The strains ND11 Malaysia and ND13 USA belonged to an independent lineage as well as distantly related to our studied strain (MSR16, MSR17) indicating that these 2 strains genome sequence might be disparate to MSR16 and MSR17 strain.
Phylogenetic relationship of identified plasmids found in the AHPND related isolates
A total of 26 V. parahaemolyticus isolates plasmid including pMSR16 and pMSR17 were selected for establishing phylogenetic relationship among the AHPND positive V. parahaemolyticus plasmid (Figure 9). The tree is unrooted as cannot determine which plasmid arose earliest. 6 plasmid including pMSR16, pVPA3-1 Vietnam, pMSR17, pVpR13-71Kb USA, pVPGX1 China, pVPE61a Thailand, were located at Cluster 1.
The phylogenetic tree showed that pMSR16 and pMSR17 were closely related with pVPA3-1. The plasmid pVPA3-1 is a Vietnam strain and its accession no is NC_025152.1. Han et al. (2015) reported that, 69 kb plasmid pVPA3-1 was identified in V. parahaemolyticus strain 13-028/A3 that can cause AHPND (9). It consists of 92 open reading frames that encode mobilization proteins, replication enzymes, transposases, virulence-associated proteins, and proteins similar to Photorhabdus insect-related (Pir) toxins. In V. parahaemolyticus, these Pir toxin-like proteins are encoded by 2 genes (pirA- and pirB-like) located within a 3.5 kb fragment flanked with inverted repeats of a transposase-coding sequence (~1 kb).
The plasmid p1937-2 China and pFORC4 South Korea: Busan strain were located at cluster 2. The plasmid pFORC14, 04-2192 located in cluster 3. The plasmid pVpR14-74Kb, pVA1, ISF-54-12, HS-22, P2, 04-2551, PMA109-5 located in the cluster 4. The plasmid S372-5, pVPUCMV, pC4602-2, RM-13-3 located at the cluster 5. The plasmid pVpR13-55Kb USA, pVpR14-56Kb USA and pVPGX2 China located at the cluster 6.
The plasmid pV110-KY498540.1 China and p1937-1-NZCP022245.1 China belonged to an independent lineage respectively. This two strains might acquire plasmids form different sources. These two plasmids are also distantly related with our studied pMSR16 and pMSR17 indicating that they might be not come from China.