To determine the prevalence and risk factors of seafood-borne V. vulnificus in Asia, a systematic review and meta-analysis were conducted following the “Preferred Reporting Items for Systematic Reviews and Meta-Analyses Protocol” (PRISMA-P) guidelines (Page et al., 2021) (Supplementary Table S1). The research question was formulated in the “population, exposure, comparator, and outcome” (PECO) format, wherein “population” indicated seafood samples and “exposure” referred to V. vulnificus. Since this systematic review and meta-analysis focused on prevalence, the term “comparator” was not applicable, hence omitted. The presence or absence of V. vulnificus in seafood samples obtained at any stage of the seafood supply chain was considered the “outcome” of each primary study.

A comprehensive literature search of four electronic databases (PubMed, Web of Science, Scopus, and Cochrane Library) was performed to identify relevant studies reporting the prevalence of V. vulnificus in seafood, published until December 31, 2022. For this purpose, an appropriate search string was developed using a combination of free text and controlled vocabulary terms such as the name of the concerned pathogen and its various pseudonyms, types of common seafood and their habitats, the names of all Asian countries, and prevalence-related keywords. The search terms were combined with Boolean operators (AND/OR), an asterisk (*), and parentheses in a meticulously selected manner to ensure the retrieval of relevant information (a detailed search string is provided in Table 1). Language restrictions were not imposed on the search, which led to the discovery of non-English articles published in Chinese, Japanese, and Korean that were translated and included in this study. All references were managed with the Endnote citation manager (Endnote X9, Clarivate, Philadelphia, PA, United States). After importing to Endnote, duplicate studies were removed. Additionally, a manual search of the reference lists in downloaded articles was conducted to identify any potential studies that might have been inadvertently overlooked during the initial search process.

After removing duplicates, the remaining articles were screened manually against predetermined inclusion and exclusion criteria (Table 1). Three independent researchers assessed the eligibility of the studies with any disagreements resolved through consensus or arbitration by an independent reviewer. An initial screening of titles and abstracts was performed, followed by a final screening of the full texts. After two-level screening, articles that fulfilled the inclusion criteria were considered eligible and selected for qualitative (systematic review) and quantitative (meta-analysis) synthesis. Non-accessible articles were obtained from the Foreign Research Information Center of Jeonbuk National University¹.

To evaluate the internal and external validity of the included studies, the risk-of-bias (RoB) assessment was carried out by three independent researchers using the critical appraisal checklist of the Joanna Briggs Institute (JBI) for prevalence studies (Martin, 2017). The checklist contains nine RoB criteria related to the target population, sampling strategy, clarity in describing study subjects, reliability and validity of diagnostic methods, data collection method, response rate, likelihood of nonresponse bias, adequacy of sample size, and appropriateness of statistical analysis. Only seven of them were evaluated because two survey-type criteria (in terms of participants’ response rates) did not apply to this study. Based on seven of the nine RoB criteria included in the checklist, each study was rated as having either a low, high, or unclear RoB. A final decision was reached through consensus among the researchers, and disagreements were resolved through arbitration by an independent reviewer. The overall RoB results of the included studies were graphically presented using the Review Manager software RevMan 5.4.1.

After validating quality of the eligible studies through RoB assessment, relevant information was collected and organized in a Microsoft Excel spreadsheet to facilitate data analysis. The extracted data included several key variables, such as author name, publication year, country of study, sample type (fish, shellfish, or other types of seafood), sampling stage (harvest or post-harvest), prevalence rate of V. vulnificus (number of total and positive samples), and detection method (biochemical tests or PCR). In cases where more than one trial was conducted in a single study (e.g., using different sample types or sampling stages), the data were extracted separately from each trial.

Given the nature of this systematic review, which included studies from various Asian regions with potentially varying conditions and populations, a random-effects meta-analysis was employed to provide a more conservative estimate of the summary effect size while accounting for the expected heterogeneity. It considers both within-study and between-study variations, making it a suitable choice for pooling data from diverse sources (Harrer et al., 2021). In the random-effects model, the pooled prevalence estimate is calculated as the weighted average of the observed prevalence estimates from individual studies. The weight assigned to each study is determined by the inverse of its variance. Consequently, the studies with smaller variances (more precise estimates) receive greater weights, while the studies with larger variances (less precise estimates) receive smaller weights. This approach ensures that more reliable estimates contribute more substantially to the overall pooled estimate, reflecting the varying levels of precision across the included studies (Lipsey and Wilson, 2001). To fulfill the assumption of normal distribution and stabilize variances, the prevalence estimates from individual studies were subjected to a logit transformation (Barendregt et al., 2013). A Generalized Linear Mixed Model (GLMM) was then employed for pooling the transformed data, and between-study variance (τ²) was calculated by the Maximum Likelihood (ML) estimator (Lin and Chu, 2020). For ease of interpretation, the prevalence rates and the corresponding 95% confidence intervals (CIs) of the logit model were back-transformed to their original forms and represented as percentages. The between-study heterogeneity was assessed using the Q test and I² statistics, which accounts for the percentage of the total variability in the effect estimates caused by true heterogeneity rather than random chance. Heterogeneity was considered significantly high if the Q test yielded a significant p-value (<0.05) and I² was greater than 50% (Higgins and Thompson, 2002; Deeks et al., 2019). The findings of the meta-analysis were visually presented as forest plots. The methodology for calculating the logit, standard error, inverse variance weight for individual studies, and the weighted average proportion using logit transformation followed the approach outlined by Wang (2023). In order to assess the robustness of the meta-analysis findings, a sensitivity analysis using the leave-one-out method was also conducted. This approach aimed to identify any outliers or influential studies that could potentially impact the pooled prevalence estimate (Wang, 2023).

Given the anticipated heterogeneity in prevalence studies, it is commonly recommended to conduct moderator analyses to delve deeper into the reasons behind between-study heterogeneity. Similarly, in this meta-analysis, six a priori-determined covariates—publication year, country, detection method, sample size, sample type, and sampling stage—were presumed to be associated with the variations in the prevalence rate of V. vulnificus across the studies. Subsequently, comprehensive moderator analyses, employing subgroup and meta-regression analyses were conducted to investigate the potential factors underlying between-study heterogeneity using the above-mentioned variables. The statistical program “R” (version 4.1.2) and accompanying R-studio (version 1.4.1106) was utilized to perform all these analyses by employing “metafor” (version 3.8–1) and “meta” (version 5.5–0) packages.

To quantify the extent of observed variability in prevalence estimates across the studies, explained by the selected variables, meta-regression analyses were conducted using the R package “MuMIn” (version 1.43.17) and “dmetar” (version 0.0.9000). Categorical variables were encoded as dummy variables to allow their inclusion in the regression models. Initially, univariate meta-regression analyses were performed for each of the six covariates independently to assess their individual impact on the between-study heterogeneity. Covariates with a p-value ≤ 0.1 were then considered for inclusion in the multivariate meta-regression model, allowing for a simultaneous examination of their impact. To enhance the robustness of multivariate meta-regression model, the multi-model inference approach was employed. This approach systematically considers and compares multiple models, each with a different combination of covariates, rather than relying on a single “best” model. Model averaging was done by assigning weight to each model based on its performance (Harrer et al., 2021). The weights assigned to each model were determined by Akaike’s information criterion (AICc). The model with the lowest AICc was identified as the best optimal model, indicating a balance between goodness of fit and the number of predictors (Field and Gillett, 2010). Subsequently, the coefficients of this model were utilized to formulate the regression equation. A plot was generated to visually represent the averaged importance of predictors across all models. This involves ranking predictors according to their relative contribution in explaining the variability across prevalence estimates. For the variables that were not included in the multiple meta-regression model, a graphical representation was utilized to visualize their impact on the prevalence estimates across the studies using the R package “ggplot2” (version 3.3.6).

Publication bias refers to the selective publication of the studies based on their statistical significance, magnitude and direction of findings (Wang, 2017). It is a threat to systematic reviews and meta-analyses, as it could lead to a biased estimate of the effect size. In this study, the likelihood of publication bias on the prevalence of seafood-borne Vibrio vulnificus in Asia was evaluated through visual inspection of the contour-enhanced funnel plot symmetry (Palmer et al., 2008). To avoid subjective inferences and enhance interpretation, publication bias was quantitatively estimated using Egger’s regression test (Egger et al., 1997). Furthermore, a trim-and-fill method was employed to generate a non-biased effect size by imputing missing studies into the funnel plot (Palmer et al., 2008).

A systematic literature search following the PRISMA guidelines was conducted to ensure transparency in reporting our research findings. Figure 1A shows the number of research articles retrieved at each stage of the process. An initial literature search of electronic databases and other sources yielded 279 documents. After removing the duplicates, 234 records remained. A further 175 studies were eliminated after the preliminary title and abstract screening, and another 21 documents were removed because of irrelevant content found during the full-text screening. Finally, 38 articles were found to be eligible for inclusion in this meta-analysis.

Before extracting data from eligible articles, RoB was assessed against the seven relevant quality criteria (C1–C7) mentioned in the JBI critical appraisal checklist to evaluate the internal and external validity of data (Martin, 2017). A descriptive justification for each judgment was presented, with figures showing the RoB assessment (Figures 1B,C). A score sheet for the quality evaluation of all included studies is also provided in Supplementary Table S2. Five of the seven examined criteria (C1, C3, C4, C6, and C7) indicated a low RoB for each study. The criteria for the appropriateness of the sample type (C1) ensured that the intended target population aligned with the inclusion requirements (Table 1). To authenticate the detection approaches, the criteria for evaluation (C3) guaranteed that each positive sample was identified efficiently, and all samples underwent the same methodology to detect the pathogen of interest (C4). Similarly, the criteria for measuring the outcome included a clear depiction of the number of positives, the number of samples tested, and/or the ability to calculate the prevalence (C6) for a specific subgroup (C7). In contrast, the RoB criteria for sampling strategy (C2) were unclear in 95% of the studies. This is primarily because random sampling was mentioned in only two studies (Amalina et al., 2019; Zangoei-Fard et al., 2020), leaving out 36 other studies. Finally, for the studies that did not report their sample size calculation, the adequacy of the used sample size was evaluated by using the formula suggested by Naing et al. (2006). By employing that approach, a minimum sample size of 150 was deemed appropriate for precise V. vulnificus prevalence estimation. Consequently, 23 of the 38 studies (61%) were judged to be at a high RoB on the adequate sample size criterion (C5) because they used less than 150 samples. Overall, all the studies fulfill at least five of the seven quality criteria, thereby warranted to be at low RoB.

The 38 cross-sectional epidemiological studies with an explicit number of samples collected for the detection of V. vulnificus prevalence were included in this meta-analysis (Supplementary Table S3). These studies were conducted between 1989 and 2022 in 11 countries located in West Asia (Turkey and Iran), South Asia (India), East Asia (China, Japan, Hong Kong, South Korea, and Taiwan), and Southeast Asia (Thailand, Malaysia, and Vietnam). Figure 2 shows a map of the Asian countries included in the meta-analysis along with the reported sample size and prevalence. Among the eligible studies, the prevalence of V. vulnificus was investigated in bivalves (oysters, mussels, clams, cockles, and scallops), crustaceans (shrimps, prawns, crabs, crayfish, and lobsters), cephalopods (squids and octopuses), fish (salmon, tuna, carp, tilapia, and sea bass), and other seafood (sea snails and sea urchins). Various stages of the seafood supply chain were categorized into two primary groups: “harvest” and “post-harvest.” Each group represents the point in the supply chain from which samples were collected. The “harvest group” encompasses the studies where seafood samples were directly obtained from natural habitats such as oceans, seas, rivers, lakes, and other water bodies. On the other hand, the “post-harvest group” comprises the studies where seafood samples were collected from seafood markets, retail shops, or other seafood stores. Along the seafood supply chain, 12 studies collected samples at the harvest stage, 23 at the post-harvest stage, and 3 studies collected samples from both stages. For the detection of V. vulnificus in seafood samples, all the included studies initiate with sample pre-enrichment (primarily utilizing alkaline peptone water), followed by pathogen cultivation and isolation on selective agar media (predominantly TCBS agar). For the species identification, 25 of the 38 studies employed PCR assays (genotypic identification) while the remaining 13 studies relied only on biochemical tests (phenotypic identification) to confirm the presence of V. vulnificus (Supplementary Table S4).

The 38 primary studies that reported the prevalence of V. vulnificus in seafood were evaluated using random effects meta-analysis. The estimated pooled prevalence was 10.47% (95% CI: 6.8–15.8), suggesting that approximately 7–16% of the seafood samples consumed in Asia contain V. vulnificus (Figure 3). Q-statistic indicated significant heterogeneity across the studies (Q = 954.48; p-value < 0.01; I² = 96%). Despite the heterogeneous prevalence rates in the primary studies, the leave-one-out analysis demonstrated that no individual study substantially influenced the outcomes (Supplementary Figure S1). The leave-one-out prevalence estimates ranged from 11% (95%CI: 8–16) to 12% (95%CI: 8–18), suggesting that the results of the meta-analysis are robust and reliable. The Subgroup analyses based on six variables—country, publication year, detection method, sample type, sample size, and sampling stage—explore the potential sources of between-study heterogeneity among the prevalence estimates (Supplementary Table S5). The subgroup analyses revealed that country, publication year, and detection method were significantly associated with variations in the prevalence estimates (p < 0.05). Regarding countries, Japan had the highest prevalence with an estimate of 47.6% (95% CI: 42.1–53.2%), followed by China (36.9, 95% CI: 23.1–53.2%), India (17.3, 95% CI: 8.3–32.5%), Thailand (16.3, 95% CI: 7.7–31.2%), and Malaysia (16.3, 95% CI: 9.7–26.3%). A prevalence rate of less than 10% was reported in Hong Kong, South Korea, Taiwan, Vietnam, Turkey, and Iran (Figure 2). Regarding publication year, the studies published before 2015 showed a higher prevalence (15.02, 95% CI: 9.02–23.97%) than those published after 2015 (6.26, 95% CI: 3.23–11.82%). In terms of detection method, the studies that included the more sensitive PCR method for the detection of V. vulnificus exhibited a higher prevalence rate (14.28, 95% CI, 8.82–22.28%) compared to those that used conventional biochemical tests (5.72, 95% CI, 2.75–11.53%). For other covariates, including sample type, sample size, and sampling stage, the results showed that although variation was present among the subgroups, the difference was not statistically significant (p > 0.05). The quantitative impact of each covariate on the between-study heterogeneity of the estimated prevalence was further assessed using meta-regression.

The relationship of between-study heterogeneity and variables of interest (country, publication year, detection method, sample type, sample size, and sampling stage) was examined using univariate and multivariate meta-regression analyses. The univariate meta-regression analysis revealed a significant relationship between country and between-study heterogeneity (R² = 18.14%, p < 0.01). The findings suggest that the country accounted for 18% of the true heterogeneity in the observed prevalence rates across the studies. Similarly, the publication year and detection method were also significantly associated with between-study heterogeneity, with R² values of 11.46 and 12.93%, respectively. Although sample size, sample type, and sampling stage contributed 9.48, 1.41, and 3.23%, respectively, to the true heterogeneity in the observed prevalence, their relationship was not significant (p > 0.05). The variables whose p-values were less than or equal to 0.1 in the univariable analyses were selected to construct a multivariable meta-regression model. Hence, the sampling stage and sample type were excluded because their p-values exceeded 0.1, whereas the other four variables were taken into account. The multivariate meta-regression model revealed that the combination of detection method, country, publication year, and sample size accounted for 48.16% of the heterogeneity across the prevalence estimates. Although this contribution was significant (p < 0.0001), the model did not entirely explain the reasons for the between-study heterogeneity. Owing to the presence of high residual heterogeneity (I² = 93.19%, p < 0.0001), it is possible that additional factors that have not been considered, may contribute to the varying prevalence estimates. The predictor importance values revealed that the detection method had the highest contribution in the regression model (r = 0.95), followed by country (r = 0.77), publication year (r = 0.72), and sample size (r = 0.45) (Figure 4A). This implies that the reasons for the variability in V. vulnificus prevalence rates in Asia are multifactorial, with the aforementioned predictors being the most important. To minimize the risk of overfitting, a model with a small number of terms and the lowest value of AICc was considered most desirable (Figure 4B). Hence, the final model was:

where logit (Y_Prevalence) denotes the logit-transformed value of prevalence and X_Country, X_{Detection method}, and X_{Publication year} are the independent study variables.

Finally, the data visualization was performed to examine the variations in the prevalence estimates across the studies based on the variables that were not considered in the multivariate meta-regression model. Given the importance of the detection method, it was also incorporated into the visualization graph to enhance the clarity of its impact on the prevalence rates among the studies. It was observed that the studies employing biochemical tests consistently reported prevalence rates below 25% over time, whereas those utilizing PCR showed a lower prevalence rate after 2015 (Figure 5A). In case of sampling stage, the subgroup analysis revealed that the prevalence of V. vulnificus in freshly harvested samples was 7.26% (CI: 3.5–14.47%), whereas a notably higher prevalence rate of 13.49% (95%CI: 8.2–21.41%) was observed in the samples collected at the post-harvest stage. Figure 5A shows that most post-harvested and all freshly harvested samples had lower prevalence rates after 2015. Regarding the sample type, the highest prevalence was detected in sea snails (27.6%), followed by cockles (19.8%), prawns (15.8%), clams (14.5%), oysters (11.7%), crayfish (11.3%), shrimps (9.4%), mussels (8%), and fish (4.6%). Conversely, V. vulnificus was not detected in octopuses or sea urchins, whereas scallops, crabs, lobsters, and squids had low detection rates (3.7, 2.2, 1.7, and 1%, respectively) (Figure 5B).

To investigate the extent of publication bias for the prevalence of seafood-borne V. vulnificus in Asia, a contour-enhanced funnel plot with different significance levels (<0.1, <0.05, <0.01) was generated with the logit-transformed proportion on the x-axis and their corresponding standard errors on the y-axis. A visual inspection of the funnel plot revealed an asymmetric distribution of the studies on both sides of the mean effect, suggesting the likelihood of publication bias (Figure 6). The presence of funnel plot asymmetry was further confirmed by Egger’s regression test that yielded a statistically significant p-value (t = −3.483, p = 0.001), suggesting the presence of publication bias for the prevalence of V. vulnificus in seafood in Asia. Applying the trim-and-fill method, 16 missing studies were identified, leading to a notable shift in the pooled prevalence estimate from 10.47% (CI: 6.8–15.8) to 28.07% (CI: 17.89–41.15). However, the shift in the prevalence estimate did not result from the inclusion of the original data, so it would not affect the validity of these findings. It rather underscores the importance of conducting additional studies to enhance the knowledge of the current prevalence status of V. vulnificus in seafood.