Impacts of Climatic Variability on Vibrio parahaemolyticus Outbreaks in Taiwan

This study aimed to investigate and quantify the relationship between climate variation and incidence of Vibrio parahaemolyticus in Taiwan. Specifically, seasonal autoregressive integrated moving average (ARIMA) models (including autoregression, seasonality, and a lag-time effect) were employed to predict the role of climatic factors (including temperature, rainfall, relative humidity, ocean temperature and ocean salinity) on the incidence of V. parahaemolyticus in Taiwan between 2000 and 2011. The results indicated that average temperature (+), ocean temperature (+), ocean salinity of 6 months ago (+), maximum daily rainfall (current (−) and one month ago (−)), and average relative humidity (current and 9 months ago (−)) had significant impacts on the incidence of V. parahaemolyticus. Our findings offer a novel view of the quantitative relationship between climate change and food poisoning by V. parahaemolyticus in Taiwan. An early warning system based on climate change information for the disease control management is required in future.


Introduction
Fish and seafood are central to the diet in Taiwan. The average per person fish consumption in Taiwan was 36.5 kg/yr in 2012, which was 2.5 times higher than the world average of 14.5 kg/yr. In Taiwan, the premier food poisoning threat over the most recent ten years has been Vibrio parahaemolyticus, and seafood is the major food type causing food poisoning (Table 1) [1,2]. V. parahaemolyticus is naturally present in the marine and estuarine environments of tropical and temperate areas. It is naturally of particular importance in the countries where seafood consumption is high [3][4][5]. V. parahaemolyticus is a Gram-negative, halophilic bacterium. There are three major manifestations of Vibrio infection: gastroenteritis, wound infection, and primary septicemia. Fatality rates were 1% for gastroenteritis, but as high as 5% for wound infections and 44% for septic disease [6]. Although the microbiological aspects of seafood safety have been studied intensively for many decades, there is still a considerable burden of food-borne illness, even in industrialized countries [7][8][9][10].
Climate change or variability is the greatest environmental challenge the world faces today [11][12][13][14]. Climate change includes changes in temperature and precipitation patterns, increased frequency and intensity of extreme weather events, ocean warming and acidification, and changes in the transport pathways of complex contaminants [15]. A substantial body of influential research indicates that climate change relates to food safety because it affects microbial ecology and growth, plant and animal physiology and host susceptibility [16][17][18]. These effects, in turn, could affect human health [11,[19][20][21][22][23][24][25].
Quantitative data have indicated that climate variation is related with infectious diseases caused by pathogenic organisms, such as Salmonella [22,26], pathogenic Escherichia coli [23], and Campylobacter [27]. Although growing numbers of researchers have considered the influences of climate change on human health caused by pathogenic organisms, very little attention has been paid specifically to V. parahamolyticus, and from the perspective of Taiwan [28]. This paper therefore serves to address this under-researched area. We aim to quantify the relationship between climate variability and V. parahaemolyticus incidence in Taiwan in order to provide a mechanism with which to establish an early warning system in the future.

Materials/Data
Monthly data for 2000-2011 of V. parahaemolyticus outbreaks (outbreak) was obtained from the Food and Drug Administration, Department of Health, Taiwan due to data access limitations. The climate data was obtained from the Central Weather Bureau, Taiwan, and were aggregated from all recording stations. They were maximum, average and minimum temperature (maxtemp, avgtemp, mintemp), maximum, average, and minimum relative humidity (maxrh, avgrh, minrh), as well as average and maximum daily rainfall (avgrf, maxrfd). Data for oceanic temperature and salinity (octemp, ocsant) were also obtained from Ministry of Science and Technology, Taiwan. A total of 144 monthly data values were obtained. The study area range is at latitude 20˝00'North to 29˝58'North and longitude 115˝00'East to 124˝99'East. Regional differences were not considered in this study because Taiwan is a small island [21].

Methods/Analysis
In order to accommodate the autocorrelation and seasonality in the time course of V. parahaemolyticus outbreak cases, we employed a Box & Jenkins [29] seasonal autoregressive integrated moving average model (ARIMA). This model captured variation, autocorrelation, long-term trends, and allowed the examination of the independent impacts of the covariates such as the climatic variables. This model has been widely used to examine the seasonality in the time course of other pathogens, such as Salmonella [22,26], pathogenic E. coli [23], and Leptospira interrogans [30]. The seasonal ARIMA model incorporates both non-seasonal and seasonal factors in a multiplicative model. One shorthand notation for the model was ARIMA (p,d,q)ˆ(P,D,Q) S . The model could be written more formally as: where y t was the outbreak variable, ε t was a white noise process, B was lag (or back) operator; non-seasonal components AR(p): , and d was the degree of non-seasonal differentiation, D was the degree of seasonal differentiation, and S was the length of the seasonal cycle.
On the left side of Equation (1) the seasonal and non-seasonal AR components multiplied each other, and on the right side of Equation (1) the seasonal and non-seasonal MA components multiplied each other. The seasonal AR components captured the seasonal pattern of the dependent variable, and non-seasonal AR components captured the lag-time structure of the dependent variable.
The outbreak variable was first tested for seasonal variation and stationarity by the Augmented Dickey-Fuller (ADF) unit root test. The orders of the ARIMA(p,d,q)ˆ(P,D,Q) S model were then determined by examining the autocorrelation function (acf) and the partial autocorrelation function (pacf). A detailed explanation of ARIMA procedures was provided in Box and Jenkins [29]. The climate variables were included in the model by calculating the cross-correlations function (ccf) to examine the correlation between the climatic variables and the V. parahaemolyticus outbreaks. All ARIMA modelling and the corresponding statistical tests were performed in the SAS/ETS 9.3 statistical software.

Results
To avoid the collinearity problem of the independent variables, we calculated the Pearson correlation coefficients between all pairs of climate variables ( Table 2). The temperature variables, maxtemp, avgtemp, and mintemp, were highly correlated (all r > 0.99) with each other, and variable avgtemp had the highest correlation with outbreaks. We therefore selected avgtemp as the temperature variable. For the same reason, maxrfd was selected as the rainfall variable and avgrh that for relative humidity. Vibrio incidence was also significantly and positively correlated with average temperature (avgtemp) (r = 0.2675), maximum daily rainfall (maxrfd) (r = 0.1402), and ocean temperature (octemp) (r = 0.1938). Average relative humidity (avgrh) and ocean salinity (ocsant) were not clearly correlated with Vibrio outbreaks.
We summarized the descriptive statistics of outbreaks and climate variables in Table 3. The average values, maximum values and minimum values of outbreak variable and climatic variables were calculated by month from 2000 to 2011. In total, 3870 V. parahaemolyticus outbreaks occurred between 2000 and 2011 in Taiwan. On average, 26.9 cases each month occurred shown in the last column of Table 3. Over 80% of incidences occurred between May and September which were the warm, humid and rainy months in Taiwan. The average temperature (avgtemp) during these months was about 3.94˝C higher than the annual average temperature, the average maximum daily rainfall was 37.36 mm more than the annual average (69.4 mm), and the average ocean temperature was 0.38˝C higher than the annual average (20˝C). In contrast, the ocean salinity during this period was 0.18 Practical Salinity Unit (psu) lower than the average (34.1 psu). Outbreaks and all these climatic variables showed annual seasonality.
Monthly incidence of Vibrio parahaemolyticus outbreaks in Taiwan from January 2000 to December 2011 compared to climatic variables for the same period is shown in Figure 1. The plot of the observed Vibrio parahaemolyticus incidence showed three major outbreaks in Taiwan (2002/4-2002/9, 2003/6-2003/9, 2008/1). The bivariate analysis between crude climatic variables and V. parahaemolyticus incidence shows that the three major outbreaks were correlated to an increase of average temperature ( Figure 1A), maximum daily rain fall ( Figure 1B), relative humidity ( Figure 1C), and a slightly increase of ocean temperature ( Figure 1D). V. parahaemolyticus incidence was not clearly correlated to monthly average ocean salinity ( Figure 1E). An annual seasonality was identified for all these climatic variables.
As a first step in ARIMA modelling, we made the response series stationary, i.e., the monthly V. parahaemolyticus outbreaks count in Taiwan. To test of stationarity, the Augmented Dickey-Fuller (ADF) unit root test was performed. The ADF test determined whether the autoregressive term had a unit root. The model, with the time trend specification, was: where, a 0 was a constant, a 2 , r the coefficients on a time trend (t) and lag 1 period of y t , P the lag order of the autoregressive (AR) process, and β i the coefficient of the lag order p of AR terms. The hypotheses of ADF test was Ho: γ = 0, H: r < 0, and MacKinnon [31] one-sided p-values was performed.
The testing results were shown in Table 4. The ADF test statistics were´9.20254,´10.8712 and´11.18945 for the models of zero mean without trend, single mean without trend and single mean with trend specification, respectively. It indicated that all of the p-values were small enough (p < 0.05) to cause to reject the null hypothesis that the series had a unit root, so that the time series of V. parahaemolyticus outbreaks was stationary.  maximum relative humidity (maxrh), minimum relative humidity (minrh), ocean temperature (octemp), and ocean salinity (ocsant); b : *, ** and *** indicate significance at the 10%, 5%, and 1% level of probability, respectively.    We further used the auto-correlation (acf) and partial auto-correlation function (pacf) plots to identify the order of the ARIMA model for the stationary series. The acf and pacf pattern of the original time-series determined the initial autoregressive (AR) and moving average (MA) order. To account for annual seasonality, we fitted several further univariate ARIMA models with different orders and subsequently excluded any models in which the residuals exhibited autocorrelation. Of which four univariate models, ARIMA(1,0,0)12, ARIMA(0,0,1)12, ARIMA(1,0,1)12 and ARIMA(1,0,0) × (1,0,0)12 (null model#1, model#2, model#3, and model#4) had showed better predictive Root Mean Square Error (RMSE) and lowest Akaike Information Criterion (AIC). The estimated models were tested for residuals autocorrelation. The presence of autocorrelation was investigated using the acf, the pacf and the correlogram reported in the Figure 3A-D. Figure 3 indicated all of the autocorrelation coefficients of the residuals did not differ significantly from zero with all acf and pacf values fell between the upper and lower bounds of 2 standard errors. All the performed diagnostic statistics indicated that the models pass all the tests. Their maximum likelihood estimation (MLE) estimated coefficients were summarized in Table 5. Of these four univariate models, ARIMA (1,0,1)12 had the best predictive RMSE and lowest AIC. We, therefore, used the model ARIMA (1,0,1)12, with AR at lag 12 and MA at lag 12, as the baseline univariate model for further comparisons. The selected model could be written as follows: We further used the auto-correlation (acf) and partial auto-correlation function (pacf) plots to identify the order of the ARIMA model for the stationary series. The acf and pacf pattern of the original time-series determined the initial autoregressive (AR) and moving average (MA) order.  We further used the auto-correlation (acf) and partial auto-correlation function (pacf) plots to identify the order of the ARIMA model for the stationary series. The acf and pacf pattern of the original time-series determined the initial autoregressive (AR) and moving average (MA) order. To account for annual seasonality, we fitted several further univariate ARIMA models with different orders and subsequently excluded any models in which the residuals exhibited autocorrelation. Of which four univariate models, ARIMA(1,0,0)12, ARIMA(0,0,1)12, ARIMA(1,0,1)12 and ARIMA(1,0,0) × (1,0,0)12 (null model#1, model#2, model#3, and model#4) had showed better predictive Root Mean Square Error (RMSE) and lowest Akaike Information Criterion (AIC). The estimated models were tested for residuals autocorrelation. The presence of autocorrelation was investigated using the acf, the pacf and the correlogram reported in the Figure 3A-D. Figure 3 indicated all of the autocorrelation coefficients of the residuals did not differ significantly from zero with all acf and pacf values fell between the upper and lower bounds of 2 standard errors. All the performed diagnostic statistics indicated that the models pass all the tests. Their maximum likelihood estimation (MLE) estimated coefficients were summarized in Table 5. Of these four univariate models, ARIMA (1,0,1)12 had the best predictive RMSE and lowest AIC. We, therefore, used the model ARIMA (1,0,1)12, with AR at lag 12 and MA at lag 12, as the baseline univariate model for further comparisons. The selected model could be written as follows: To account for annual seasonality, we fitted several further univariate ARIMA models with different orders and subsequently excluded any models in which the residuals exhibited autocorrelation. Of which four univariate models, ARIMA(1,0,0) 12 , ARIMA(0,0,1) 12 , ARIMA(1,0,1) 12 and ARIMA(1,0,0)ˆ(1,0,0) 12 (null model#1, model#2, model#3, and model#4) had showed better predictive Root Mean Square Error (RMSE) and lowest Akaike Information Criterion (AIC). The estimated models were tested for residuals autocorrelation. The presence of autocorrelation was investigated using the acf, the pacf and the correlogram reported in the Figure 3A-D. Figure 3 indicated all of the autocorrelation coefficients of the residuals did not differ significantly from zero with all acf and pacf values fell between the upper and lower bounds of 2 standard errors. All the performed diagnostic statistics indicated that the models pass all the tests. Their maximum likelihood estimation (MLE) estimated coefficients were summarized in Table 5. Of these four univariate models, ARIMA (1,0,1) 12 had the best predictive RMSE and lowest AIC. We, therefore, used the model ARIMA (1,0,1) 12 , with AR at lag 12 and MA at lag 12, as the baseline univariate model for further comparisons. The selected model could be written as follows: where, y t represented the V. parahaemolyticus outbreaks. In order to incorporate the climate variables as input variables to the baseline univariate model, we first used ccf to examine the correlations between V. parahaemolyticus outbreaks and the environmental series. Figure 4 indicated that there were significant correlations (based on the two standard error limits) of outbreaks with monthly average temperatures at lag 0 to 8, monthly average relative humidity at lag 8 to 10, monthly maximum daily rainfalls at lag 7 and ocean temperature at lag 0. It was therefore evident that the original time-series was correlated with these climatic factors. Subsequently, the baseline univariate model (ARIMA (1,0,1) 12 ) with one or more climatic variables included were estimated. Several models were computed and only models with statistically significant coefficient were selected, ensuring the non-autocorrelation of residuals at 5% significant level. The estimated coefficients of ARIMA(1,0,1) 12 for V. parahaemolyticus outbreak series with climatic variables were shown in Table 6. For these covariate models, the best fit was obtained from ARIMA(1,0,1) 12 with average temperature, current and one month previous maximum daily rainfall, current and 9 months previous relative humidity, ocean temperature, and ocean salinity of 6 months previous as covariates had the best prediction RMSE and lowest AIC. The selected model could be written as follows: where, y t represented the V. parahaemolyticus outbreaks, y t´12 12 months previous outbreak, avgtemp t current average temperature, maxrfd t current maximum daily rain fall, maxrfd t´1 one month previous maximum daily rain fall, avgrh t current average relative humidity, avgrh t´9 9 months previous average relative humidity, octemp t current ocean temperature, and ocsant t´6 ocean salinity of 6 months previous. Comparing the model with the baseline univariate model (ARIMA(1,0,1) 12 ), we found that including the environmental input series improved the AIC by 1.57% and the predictive RMSE by 2.54% compared to the baseline model.
where, yt represented the V. parahaemolyticus outbreaks. In order to incorporate the climate variables as input variables to the baseline univariate model, we first used ccf to examine the correlations between V. parahaemolyticus outbreaks and the environmental series. Figure 4 indicated that there were significant correlations (based on the two standard error limits) of outbreaks with monthly average temperatures at lag 0 to 8, monthly average relative humidity at lag 8 to 10, monthly maximum daily rainfalls at lag 7 and ocean temperature at lag 0. It was therefore evident that the original time-series was correlated with these climatic factors. Subsequently, the baseline univariate model (ARIMA (1,0,1)12) with one or more climatic variables included were estimated. Several models were computed and only models with statistically significant coefficient were selected, ensuring the non-autocorrelation of residuals at 5% significant level. The estimated coefficients of ARIMA(1,0,1)12 for V. parahaemolyticus outbreak series with climatic variables were shown in Table 6. For these covariate models, the best fit was obtained from ARIMA(1,0,1)12 with average temperature, current and one month previous maximum daily rainfall, current and 9 months previous relative humidity, ocean temperature, and ocean salinity of 6 months previous as covariates had the best prediction RMSE and lowest AIC. The selected model could be written as follows: where, yt represented the V. parahaemolyticus outbreaks, yt-12 12 months previous outbreak, avgtempt current average temperature, maxrfdt current maximum daily rain fall, maxrfdt-1 one month previous maximum daily rain fall, avgrht current average relative humidity, avgrht-9 9 months previous average relative humidity, octempt current ocean temperature, and ocsantt-6 ocean salinity of 6 months previous. Comparing the model with the baseline univariate model (ARIMA(1,0,1)12), we found that including the environmental input series improved the AIC by 1.57% and the predictive RMSE by 2.54% compared to the baseline model.        a The four univariate models can be represented as follows: model#1, ARIMA(1,0,0) 12 : (1´Φ 1 B 12 )y t = u + ε t , model #2, ARIMA(0,0,1) 12 : y t = u + (1 + ΘB 12 ) ε t , model #3, ARIMA(1,0,1) 12 : (1´Φ 1 B 12 )y t = u + (1 + ΘB 12 )ε t , model #4, ARIMA(1,0,0)ˆ(1,0,1) 12 : (1´Φ 1 B 12 )(1´φB)y t = u + (1 + ΘB 12 )ε t . where, y t is outbreak; ε t is a white noise process; u is intercept; Φ, Θ and φ are estimate parameters, and B is lag (or back) operator; b *** indicate significance at the 1% level of probability, respectively. The best fit model indicated that average temperature, ocean salinity of 6 months previous, and ocean temperature were all significantly and positively related to V. parahaemolyticus outbreaks. Thus outbreaks might increase when increasing in current temperature or ocean temperature or, when increasing in ocean salinity of 6 months previous. However, there were significant negative relationships of outbreaks with maximum daily rainfall of a month previous or relative humidity of 9 months previous. This relationship meant that increased heavy rain or increased relative humidity 9 months previously may reduce current outbreaks.

Key Findings
The results from our models were consistent with expectations and with results in the literature. The variables associated with food-borne diseases were ambient temperature, ocean temperature, ocean salinity, and rainfall. Ambient temperature, ocean temperature, and ocean salinity were positively related to V. parahemolyticus outbreaks. Rainfall was negatively related to V. parahemolyticus outbreaks.
Theoretically, temperature was likely to influence the growth of V. parahaemolyticus. V. parahaemolyticus can grow well at below 10˝C, moreover, its generation time can be as fast as 10 min at~37˝C [32]. Our empirical evidence showed that environmental temperature was the major factor determining the seasonality of growth and the geographical distribution of V. parahaemolyticus in shellfish [33,34]. The results from our models were largely consistent with expectations and results published in the literature. For example, there was a significant positive correlation between the mean temperature of the previous month and the number of salmonellosis notifications in the current month in five Australian cities. The increase in notifications was 4%-10% for every degree of temperature rise. Tam et al. [27] also investigated the relationship between ambient temperature and Campylobacter enteritis using time-series analysis to study short-term associations between temperature and number of Campylobacter reports in England and Wales. They also adjusted for long-term trend and seasonal patterns. They found a linear relationship between mean weekly temperature and reported C. enteritis.
In addition, our results also confirmed the importance of salinity as a climate predictor of seafood poisoning outbreaks. Compared to V. vulnificus, V. parahaemolyticus could tolerate higher salinity levels [35,36]. V. parahamolyticus could survive at salinity ranging from 10 to 34 ppt, with 23 ppt as its optimum salinity [28]. Although recent researches had considered impacts of climate change on salinity in costal and marine ecosystems [15,37], very little empirical evidence had been given on influences of salinity on seafood poisoning outbreaks. Recently, Young et al. cconducted meta-analysis to exam impacts of salinity on V. parahaemolyticus in oysters at harvest and in harvest waters [38]. They found no consistent relationship for water salinity. They further explained this might be due to poor reporting of study sampling methods and quantitative outcome data from collected articles. We thus encourage further research to investigate the impact of salinity in other regions or in other waterborne pathogens.
Rainfall was negatively related to V. parahemolyticus outbreaks as in covariate #3 and covariate #4. Deepanjali, Kumar and Karunasagar [34] also found that V. parahaemolyticus abundance was high during the dry season and low after rainy seasons. This might result from pathogen concentrations in estuaries being diluted by heavy rain thus reducing contamination of oysters and other filter feeding shell fish. In both Brisbane and Townsville rainfall was also a significant factor in regression models of climate effects on cases of Salmonella infection [26].
Furthermore, our findings also indicated that ambient temperature, ocean temperature, ocean salinity had higher impacts than rainfall on microbial seafood safety. This echoed reports in the literature showing that temperature change plays an important role in food-borne diseases [22,23,26,27]. In distinction to other studies of food pathogen outbreaks, we also proved here the importance to V. parahaemolyticus outbreaks not only of ambient temperature, but also of ocean factors such as ocean temperature and ocean salinity. Marques et al. [11] contend that a combination of environmental and genetic factors might play a key role in the presence or absence of virulent Vibrio spp. because it may alter the deposition of material in bivalve shells or the composition of plankton exoskeletons. This will in turn directly change the habitat of Vibrio spp., so forcing the organism to adapt genetically and produce strains with great virulence. We likewise demonstrated that knowing the effects of multiple climate factors on food-borne diseases was necessary for an intensive examination of accuracy in predicting disease. Moreover, effort should also be devoted to individualize more precisely the consequences of these interactions on V. parahaemolyticus infections in order to better identify control and mitigation measures [39].
Overall, increases of seafood contamination by V. parahaemolyticus in Taiwan due to climate changes could be explained as below. V. parahaemolyticus is part of the natural flora of estuarine and coastal marine environments. Climate change will likely influence the vulnerability of estuaries to eutrophication in several ways, including changes in temperature, sea level, and exchange with the coastal ocean and salinity [15,40,41]. This then will influence aquatic animals, which are vulnerable to climate change because their related metabolic processes are influenced by water temperature, salinity, and oxygen levels [42]. This in turn may favour a group of potentially emerging microbiological pathogens, the marine Vibrios, which are a genus of thermodependent bacteria which thrive in naturally warm sea water. Furthermore recent studies regarding climate impacts on marine systems in Taiwan has provided some evidence. For example, Chang et al. [43] and Lee et al. [44] observed that extreme weather and marine environmental changes induced by climate change could harm the marine fish population and aquaculture. Lu et al. [45] also observed that increased sea surface temperatures could causes fluctuations in the presence of cold-water and warm-water fishes and in the timing of fishing seasons in coastal zones of the Kuroshio Current and China Coastal Current. Chou et al. [21] also proved that climatic variations could influence diarrhea-associated morbidity in Taiwan. Thus, we can speculate that the increased seafood contamination by V. parahaemolyticus in Taiwan could be caused by climate changes.

Limitations
This research made an important contribution to the academic literature and provided the potential for positively influencing risk management practice; there were nonetheless limitations that provide opportunities for further research. First, our results were only for Taiwan. To generalise our findings we encouraged the implementation of similar studies in other countries. Second, we evaluated only temperature, humidity, and rainfall as predictors. We support, therefore, the inclusion in future studies of other oceanic factors, such as ocean turbidity, dissolved oxygen, and pH. For example, Parveen et al. [33] indicated that dissolved oxygen might increase V. parahaemolyticus abundance. Third, our analysis found that key predictors were related to the number of cases of V. parahaemolyticus infection. We suggested therefore that future studies should investigate the individual effects of these predictors on fishery production, fishery manufacturing, and the distribution system [10,46].

Conclusions
This research provided empirical evidence for the relationship between V. parahaemolyticus outbreaks and climatic change factors. Our results showed that average temperature, ocean salinity of 6 months previous, and ocean temperature were all significantly and positively related to V. parahaemolyticus outbreaks. However, there were significant negative relationships of outbreaks with maximum daily rainfall of a month previous or relative humidity of 9 months previous. Our findings also indicated that ambient temperature, ocean temperature, ocean salinity had higher impacts than rainfall on microbial seafood safety. Overall, our findings predicted that, in Taiwan, food poisoning caused by V. parahaemolyticus was related with climate factors. It is hoped that the findings of this research will help to guide public health and community interventions to protect Taiwanese consumers. Thus, future research into the health impacts of extreme weather events, and early warning and response tools is required.