Introduction

The Gram-negative, rod-shaped, halophilic bacterium Vibrio parahaemolyticus is usually found in marine and estuarine environments (Ellingsen et al. 2008). It is one of the main bacterial pathogens responsible for food-borne outbreaks and acute gastroenteritis. People are typically infected by eating raw or undercooked seafood (Pal and Das 2010). Some virulence factors may be involved in the pathogenesis, such as thermostable direct hemolysin (TDH), thermostable direct hemolysin-related hemolysin (TRH), type 3 secretion systems 1 2 (Makino et al. 2003).

Vibrio parahaemolyticus has been identified as the major cause of seafood-associated diarrhea and gastroenteritis since the pandemic O3:K6 serotype first emerged in Asia in 1996 (Okuda et al. 1997). Since then, the O3:K6 serotype and its relevant serovariants have been reported in many countries, including the United States, Mexico, Mozambique, Spain, Russia, China and elsewhere (Nair et al. 2007).

Currently, V. parahaemolyticus is the predominant cause of bacterial food poisoning in China (Wang et al. 2007), exerting a huge burden on national health care. Therefore, it is necessary to explore its molecular characteristics and epidemiology for effective prevention and control of V. parahaemolyticus. In recent years, a number of relevant studies on V. parahaemolyticus have been conducted in many Chinese cities, such as Beijing, Jiangsu, Zhejiang, Fujian, Guangdong and Taiwan (Li et al. 2014; Chao et al. 2009; Yu et al. 2011; Fan et al. 2013; Han et al. 2012). However Weihai, a typical coastal city with rapid development in north China, either lacks such data or has no data at all.

With the development of molecular biology, various approaches have been applied to characterise V. parahaemolyticus (Foxman et al. 2005), including pulsed-field gel electrophoresis, arbitrarily primed polymerase chain reaction, microarray-based comparative genome hybridization and multi-locus sequence typing (MLST) (González-Escalona et al. 2008; Wong et al. 2007; Broberg et al. 2011).

In the present study, we selected the MLST scheme developed by González-Escalona et al. (2008). It has a better accuracy in tracking genetic differences, population genetics and phylogenetic relatedness of different strains (Banerjee et al. 2014; Ellis et al. 2012). MLST has the advantages of simple operation, time-saving, a public database and easy comparison between different laboratories. It also contributes to surveys of infection sources and transmission routes (Pérez-Losada et al. 2013).

Taking full advantage of the PubMLST database, we sought to elucidate the prevalence and genetic diversity of V. parahaemolyticus circulating in Weihai at the nucleotide and peptide levels. In summary, our data not only increased the richness of the database, but also identified difference from other provinces.

Materials and methods

Sampling of V. parahaemolyticus strains

A total of 40 V. parahaemolyticus isolates (including 25 clinical isolates and 15 environmental isolates) were collected for the MLST analysis. Among them, 28 strains were isolated from different regions in Weihai, China. The other 12 isolates were obtained to enrich the database’s diversity, including three from Japan, one from the United States and eight from Thailand; these strains were collected and preserved by the General Hospital of the Jinan Military Region of the People’s Liberation Army. All isolates were identified by VITEK2 system and 16S rRNA gene sequencing. Detailed information including isolation year, site and source are listed in Table 1. The isolates were stored in Luria–Bertani broth containing 30% glycerol at − 80 °C.

Table 1 Sources of the isolates analysed

PCR amplification and sequencing

The PCR amplification protocols and relevant primer information are available at http://pubmlst.org/vparahaemolyticus/info/protocol.shtml. All isolates were PCR amplified and sequencing conducted by Thermo Fisher Scientific using M13 primers.

Diversity of pSTs

Seven nucleotide sequences were queried against the V. parahaemolyticus PubMLST database to determine the allelic profile and ST of each isolate. When a new allelic profile or ST number was discovered, the allelic nucleotide sequences were emailed to http://www.narjol.gonzalez-escalona@fda.hhs.gov, and then they were assigned a novel numerical identifier for allelic profile or ST.

As described by others (Urmersbach et al. 2014), each locus peptide sequence was queried against PubMLST database to be assigned to an allelic profile and the seven different allelic profiles led to a pST at the peptide level. The loci were labeled with the prefix ‘p_’, and scheme was called AA-MLST, which allows an analysis at the phenotypic level.

Isolate descriptions (such as place, source and isolation year), allelic profiles and sequencing traces (unique alleles) have been uploaded and made available in the PubMLST database.

Locus statistics

DnaSP Version5 (Librado and Rozas 2009) was implemented to calculate the number of alleles, number of polymorphic sites and nucleotide diversity for each locus. START V2 (Jolley et al. 2001) was implemented to calculate dN/dS through the Nei and Gojobori method. The standardised index of association (I SA ) was also caculated by START V2.

Assignment to CCs

The CCs were assigned by eBURST version3 using STs as above-described. As reported previously (Turner et al. 2013), singletons were defined when one ST differed at two or more alleles from every other ST (Feil et al. 2004); a single locus variant (SLV) was named when two STs differed by only one single locus; and a double locus variant (DLV) was defined when two STs differed by two loci (González-Escalona et al. 2008). Likewise, pSTs were assigned to CCs by goeBURST algorithm as previously described (Urmersbach et al. 2014).

The “population snapshot” analysis was also implemented with the group definition zero of seven loci shared and to create a full-(AA)MLST, where all (p)STs were connected (Francisco et al. 2009).

Phylogenetic analysis

To estimate the genetic distance of the analysed isolates, ME trees for the in-frame concatenated sequences (dnaE-gyrB-recA-dtdS-pntA-pyrC-tnaA) of each (p)ST were constructed with Mega 6 software using the Kimura two-parameter model, and the number of bootstrap replications was 1000 (Han et al. 2014).

Detection of virulent determinants and antimicrobial susceptibility tests

PCR for the detection of pathogenicity-related genes (TDH and TRH) was performed on all strains independently, using the primers described by Khouadja et al. (2014). Antimicrobial susceptibility was tested using the VITEK2 system.

Comparison of clinical isolates between Weihai and other regions in China

Using the V. parahaemolyticus PubMLST database, we searched the clinical isolates from China, and some isolates with incomplete information were excluded. Subsequently, we attempted to find the similarities and differences in clinical isolates between Weihai and other regions in China.

Results and discussion

A total of 40 isolates were collected in the present study, including 25 clinical isolates and 15 environmental isolates. Among them, 28 were isolated from different regions in Weihai, China. The other 12 isolates were isolated in other countries, including three from Japan, one from the United States and eight from Thailand (Table 1).

Diversity of peptide STs (pSTs)

Table 2 summarises data on allelic profiles and pSTs of the 40 isolates at the nucleotide and peptide levels. 23 STs were identified by MLST analysis. Among these STs, 18 STs (20 isolates) and one allele (of tnaA) were identified for the first time, including seven (29.2%) clinical isolates and 13 (54.2%) environmental isolates. New STs were mostly identified based on the combination of alleles. Most STs were represented by a single isolate, while just four STs ranged from two to 10 isolates. The most frequently identified STs were ST3 (25%) and ST332 (17.5%) in clinical isolates, whereas no predominant ST was observed from environmental isolates.

Table 2 The allelic profiles and (p)STs of the analysed isolates

After the in-frame nucleotide sequences were translated into peptide sequences, 13 distinct pSTs were acquired with a frequency ranging from 2.5 to 37.5%. This allows an analysis at the phenotypic level, and can be considered when non-synonymous substitutions of nucleotides lead to a different amino acid. Among these pSTs, there were five (38.5%) new pSTs, all belongs to environmental isolates. pST2 (37.5%, 15/40) and pST3 (25%, 10/40) were the two predominant pSTs. As reported previously (Han D et al. 2014), one particular pST can contain a large amount of nucleotide STs. In this study, we found that pST2 contained five different STs (ST3, ST479, ST1533, ST1540 and ST1547) and pST3 had four STs (ST332, ST1532, ST1534 and ST1546).

Locus statistics

Table 2 summarises the data on diversity of the seven loci, including fragment size, number of alleles, number of polymorphic sites, nucleotide diversity (π) and dN/dS ratio. The number of alleles observed for each ST ranged from 12 (tnaA) to 17 (gyrB), and the most frequently detected alleles at each locus were recA 19 and 141 (10), gyrB 4 (10), dnaE 3 (13), dtdS 4 (10), pntA 29 (11), pyrC 4 (10) and tnaA 22 (10), respectively. In addition, the number of alleles for each pST ranged from 1 (tnaA/recA) to 5 (pyrC). One allele was dominant for most alleles (~ 80%), with the exception of p_dnaE (accounting for 50.0%) and p_pyrC (accounting for 62.5%). The number of polymorphic nucleotide sites ranged from 83 (pyrC) to 159 (gyrB), while the number of polymorphic peptide sites ranged from 2 (dtdS) to 24 (gyrB). The gyrB allele possessed the largest number of polymorphic sites at both the nucleotide and peptide levels.

The nucleotide diversity (π) ranged from 0.01909 (tnaA) to 0.04107 (dnaE). The dN/dS ratio was close to zero for all the seven loci (Table 3), suggesting that nonsynonymous sites were evolving more slowly than synonymous sites, and the population was mainly affected by negative selection (purifying selection) during the evolutionary process (Jolley et al. 2001; González-Escalona et al. 2008; Han et al. 2015). This finding could explain the phenomenon that the number of alleles was reduced and the loci were primarily dominated by a single allele when analysed at the peptide level (Pérez-Losada et al. 2006; Turner et al. 2013). For example, 12 different tnaA alleles at the nucleotide level were transformed into one p_tnaA at the peptide level, which was reflected by a low dN/dS value of 0.0073. In this study, the calculated I SA value was 0.762 (P < 0.01), indicating that these alleles were in linkage disequilibrium or not randomly distributed in the V. parahaemolyticus population (Lüdeke et al. 2015). To some extent, our findings supported the hypothesis that the population structure of V. parahaemolyticus followed the epidemic model of clonal expansion (Yu et al. 2011).

Table 3 Diversity of the seven loci in the 40 V. parahaemolyticus isolates

Assignment to clonal complexes (CCs)

The global ST database was referenced to query our profiles and the results showed that the 23 STs belonged to four CCs, six groups and nine singletons (ST1536, ST1537, ST1539, ST1541, ST1542, ST1543, ST1545, ST1546 and ST1550) (Fig. 1). The CCs included CC3 (ST3), CC83 (ST1 and ST1547), CC332 (ST332, ST1532 and ST1534), CC345 (ST1533) and CC1256 (ST1256). In our study, CC3 was the most populated CC (10 isolates), followed by CC332. Among these CCs, this analysis revealed four SLVs (ST332 and 1532; ST332 and 1534; ST1540 and 1548; ST1 and 1547) and one DLV (ST1532 and ST1534). The nine singletons demonstrated that most of the STs identified in our study shared less than six alleles.

Fig. 1
figure 1

The CCs assigned by eBURST version 3 using STs

However, neither the relationships between the nine singletons themselves nor the relationships between the singletons could be deduced here. This suggests that goeBURST had limited application for identifying related isolates at the nucleotide level (Han et al. 2014).

To challenge this, a “population snapshot” analysis was implemented at the peptide level (Fig. 2). The result showed that all pSTs (except for pST10, pST241 and pST240) led to a single complex founded by pST2 (founder) and pST1 (subgroup founder). Therefore, the relationship among the isolates became more reliable when analysed at the peptide level compared with the nucleotide level.

Fig. 2
figure 2

The CCs assigned by eBURST version 3 using pSTs

Phylogenetic analysis

A minimum evolution (ME) tree was constructed using the concatenated nucleotide sequences (Fig. 3). From the ME tree, we can see that the isolates were grouped into two clusters. Cluster II only contained four environmental isolates from Weihai. Cluster I consisted of most isolates originating from different areas and sources, showing a high degree of genetic diversity. Isolates belonging to the same CCs (CC332, CC345, CC83, CC3) or groups by eBURST were also close to each other in the ME tree. The nine singletons were spread over the branches. From the ME tree, Env-WH4 and Env-WH5 clustered together (100% bootstrap support), while the eBURST analysis identified them both as singletons. There was also a high similarity between Cli-JP1 and Cli-US1 strains, which were collected from different countries. Generally speaking, results between the eBURST and ME tree analyses were consistent. However, the ME tree provided an improved phylogenetic analysis and elucidated some genetic relationships that were not observed by eBURST.

Fig. 3
figure 3

The ME tree constructed from analysis of concatenated nucleotide sequences

The ME tree was also implemented based on the peptide sequences. Data showed that 13 pSTs were grouped into one major cluster, except for Env-WH3 (Fig. 4). The evolutionary distance was 0.017930.

Fig. 4
figure 4

The ME tree constructed from analysis of concatenated peptide sequences

Detection of virulent determinants and antimicrobial susceptibility tests

20% (8/40) V. parahaemolyticus strains harboured both the TDH and TRH genes, including 7 clinical isolates and 1 environmental isolate. 52.5% (21/40) strains only contained the TDH gene, while 27.5% (11/40) strains harboured neither of these genes. The antimicrobial susceptibility tests showed all strains were resistant to ampicillin and sensitive to ampicillin/sulbactam, piperacillin/tazobactam, cephalosporins and quinolone antibiotics.

Comparison of clinical isolates between Weihai and other regions in China

Clinical isolates from China were collected via the MLST database (http://pubmlst.org/vparahaemolyticus/). Next, we screened out them for complete sequence information. Finally, 248 strains were retained and queried against 157 STs. These isolates were collected from Taiwan, Hongkong, Zhejiang, Shanghai, Beijing and elsewhere. Among them, ST3 was predominant and found in 53 isolates (21.4%), followed by ST1464 (8, 3.2%) and ST69 (5, 2.0%). In the Weihai area, we collected a total of 21 clinical isolates. Similarly, the largest number of isolates (10, 47.6%) was of ST3, followed by ST332 (7, 33.3%).

Conclusions

Weihai is a coastal city located in the east of Shandong Peninsula, facing to the Yellow Sea from three sides. Because of its geographic position and seafood-eating habits, food-borne outbreaks and acute gastroenteritis have often occurred via eating raw or undercooked seafood polluted by V. parahaemolyticus. Therefore, it’s a necessity to understand the epidemiology and genetic diversity of V. parahaemolyticus strains in this region.

An MLST analysis based on a set of 10 loci, showing a high level of nucleotide diversity and a detailed analysis of recombination, has been reported (Yan et al. 2011). In the present study, we employed a protocol reported by González-Escalona et al. (2008), which has a public database and facilitates comparisons between different laboratories (Lüdeke et al. 2015). Moreover, we identified 23 STs and found 18 novel STs, which substantially contributed to the diversity and richness of the MLST database. The environmental isolates were entirely new found STs, mostly due to new combinations of alleles. Previous studies have confirmed that V. parahaemolyticus populations in a restricted locality can be as diverse as the entire worldwide pubMLST database (Johnson et al. 2009). Our results supported this observation. However, the diversity was decreased when an analysis based was carried out on the peptide sequences, as only 13 pSTs were generated. In addition, pST2 and pST3 were the most prevalent and most were found in clinical isolates, which is in agreement with previous reports (Urmersbach et al. 2014).

The eBURST and ME tree analyses were largely consistent with each other. The isolates belonging to four CCs and groups were also clustered together in the ME tree. Some studies have reported that the isolates belonging to the same CCs may be divided into different clusters (Ansaruzzaman et al. 2005). Such discrepancy could be attributed to different analytical approaches. The goeBRUST algorithm is an allelic profile-based approach, and only SLVs are assigned to a single CC or group. Therefore, it could quickly find CCs and their ancestral STs. However, the ME tree is established based on sequences, and sequences with fewer differences are clustered together, providing an improved resolution for revealing evolutionary relationships regardless of CCs, groups or singletons. Therefore, a more reliable conclusion could be effectively acquired by combination of these two methods.

We also compared the isolates from Weihai with those from other regions of China. As reported globally, ST3 is largely responsible for V. parahaemolyticus outbreaks (Lüdeke et al. 2015) and is the main causative agent for the pandemic CC in most regions of China, including Weihai. 74% and 94% of V. parahaemolyticus infections were caused by ST3 in Beijing (Fan et al. 2013) and Zhejiang (Han et al. 2012), respectively, which are higher than that found in Weihai (47.6%), this perhaps reflects the small sample size studied here. In Weihai, the second most frequent ST was ST332. However, Han et al. have discovered that the second most frequent allelic profile includes ST305, ST8, ST8, ST301, ST670, ST432, ST91, ST1034, ST69 and ST265 in Guangdong, Jiangsu, Taiwan, Shanghai, Zhejiang, Beijing, HongKong, Fujian and Shenzhen (2014), respectively. Based on the related literature, it seems that ST332 may be restricted to the eastern coast of Shandong. But, as Han et al. (2014) reported, the V. parahaemolyticus database does not contain all the discovered isolates and it only contains reported from voluntary uploading of laboratory data. This may have some influence on acquisition of truly representative data. In Weihai, we, for the first time, found four isolates (ST1532, ST1533 and ST1534), which were associated with V. parahaemolyticus gastroenteritis. Interestingly, Turner et al. have discovered that some environmental isolates belong to ST3 by analysing a set of strains from the United States (Turner et al. 2013). Likewise, we also obtained two environmental isolates belonging to ST1 and ST 479. They were also found from clinical sources in the public database. This suggests that pathogenic strains may emerge from the environmental condition, showing a higher potential virulence to people.

In conclusion, based on MLST analysis, ST3 and ST332 were the most epidemic clones, which accounted for a large proportion of the pathogenic strains in the coastal city of Weihai. ST332 might be a region-specific ST, which needs to be confirmed by further analysis. Some newly identified STs indicated that the genes underwent recombination frequently (Yan et al. 2011). However, the study still needs some improvement. For example, we only collected 40 isolates and this small sample size limits to some extent the conclusions that can be drawn from the findings. The MLST technique employed is also unsuited for detecting significant differences in genome inversions, transposons and plasmids. Application of pulsed-field gel electrophoresis profiling may remedy this, so further studies are needed.