Comprehensive Dynamic Influence of Multiple Meteorological Factors on the Detection Rate of Bacterial Foodborne Diseases under Spatio-Temporal Heterogeneity

Foodborne diseases are a critical public health problem worldwide and significantly impact human health, economic losses, and social dynamics. Understanding the dynamic relationship between the detection rate of bacterial foodborne diseases and a variety of meteorological factors is crucial for predicting outbreaks of bacterial foodborne diseases. This study analyzed the spatio-temporal patterns of vibriosis in Zhejiang Province from 2014 to 2018 at regional and weekly scales, investigating the dynamic effects of various meteorological factors. Vibriosis had a significant temporal and spatial pattern of aggregation, and a high incidence period occurred in the summer seasons from June to August. The detection rate of Vibrio parahaemolyticus in foodborne diseases was relatively high in the eastern coastal areas and northwestern Zhejiang Plain. Meteorological factors had lagging effects on the detection rate of V. parahaemolyticus (3 weeks for temperature, 8 weeks for relative humidity, 8 weeks for precipitation, and 2 weeks for sunlight hours), and the lag period varied in different spatial agglomeration regions. Therefore, disease control departments should launch vibriosis prevention and response programs that are two to eight weeks in advance of the current climate characteristics at different spatio-temporal clustering regions.


Introduction
Foodborne diseases are illnesses that result from eating food that has been contaminated with bacteria or other pathogens such as viruses or parasites [1]. Bacterial infection is the most common type of foodborne disease, especially in the summer. Foodborne diseases are a serious global health issue and cause substantial medical costs and productivity losses worldwide [2][3][4][5]. According to data collected by the Foodborne Disease Outbreaks Surveillance System during the years 2003-2017, there were 235,754 foodborne illnesses and 1457 deaths reported to the Centers for Disease Control and Prevention (CDC), China [6]. A total of 11.3% of the 13,307 outbreaks with known etiologies were caused by Vibrio parahaemolyticus. This ranked second among foodborne disease outbreaks after poisonous mushrooms. All types of food are potentially contaminated with bacteria. Contamination can occur during food production, harvesting, processing, storage, shipping, and preparation [7]. The source of contamination varies, and most bacteria that cause food-poisoning prefer hot and humid environments [8]. Therefore, it is meaningful to study 2. Methodology

Study Area and Data Sources
The Zhejiang Province is situated on the southeast coast of China and the southern wing of the Yangtze River Delta (Figure 1). Zhejiang is a major marine and fishery prov ince in the country with a long history of fishing. It is located in a subtropical monsoon climate area with different types of topographies across the entire region. The high tem peratures and humidity in summer and autumn are suitable for the growth and repro duction of pathogenic bacteria. The incidence of foodborne diseases in the province has increased in recent years. V. parahaemolyticus was the first pathogen to cause microbia foodborne disease outbreaks in Zhejiang and accounted for 58.41% of the total bacteria outbreaks [33]. The database of bacterial foodborne diseases caused by V. parahaemolyticus from January 2014 to December 2018 recorded the date of onset, sex, occupation, and cur rent address. The dataset source was the Foodborne Disease Surveillance And Reporting System of Zhejiang Province, which collects case reports from 101 sentinel hospitals cov ering 89 county-level jurisdictions in the province. There is a high incidence of bacteria foodborne diseases in this region. However, there are few studies analyzing the dynamic relationship between bacterial foodborne diseases and meteorological factors in Zhejiang Province. Meteorological data were provided by the China Meteorological Data Network, in cluding county precipitation, average wind speed, average temperature, average relative humidity, daily minimum temperature, daily maximum temperature, and sunlight hours In this study, the daily meteorological value of each station was interpolated by inverse distance weighted interpolation (IDW) and then fitted with the centroid of each county. Meteorological data were provided by the China Meteorological Data Network, including county precipitation, average wind speed, average temperature, average relative humidity, daily minimum temperature, daily maximum temperature, and sunlight hours. In this study, the daily meteorological value of each station was interpolated by inverse distance weighted interpolation (IDW) and then fitted with the centroid of each county.

Methodological Framework
The methodological framework is illustrated in Figure 2. The spatio-temporal pattern of bacterial foodborne diseases was initially identified using a multivariate time series analysis and spatio-temporal scanning statistics. Then, the influence of high-order multicollinearity was reduced using a PCA analysis, and the effects of seven meteorological factors on bacterial foodborne diseases were explored. Finally, the VAR model was established to determine the relationship between climate variables, and the time series of bacterial foodborne diseases was examined in the different regions identified in the first step.

Methodological Framework
The methodological framework is illustrated in Figure 2. The spatio-temporal pattern of bacterial foodborne diseases was initially identified using a multivariate time series analysis and spatio-temporal scanning statistics. Then, the influence of high-order multicollinearity was reduced using a PCA analysis, and the effects of seven meteorological factors on bacterial foodborne diseases were explored. Finally, the VAR model was established to determine the relationship between climate variables, and the time series of bacterial foodborne diseases was examined in the different regions identified in the first step.

Multivariate Time Series Analysis
A time series analysis describes the changing characteristics of time series in the past and predicts the trends of the future time series; this is an important method for studying the dynamic characteristics, periodic characteristics, and correlation of statistical indicators. Scholars usually use historical data on the incidence rate of vibriosis to predict future incidence rates. However, the incidence rate is related to historical values and also depends on the values of other relevant variables in the same time series, such as temperature, relative humidity, and other climate variables. The covariance and correlation coefficients can be used to describe the degree of correlation between the two time series variables. A correlation analysis determines the interdependence between two or more variables to measure the degree and direction of the correlation between variables and the

Multivariate Time Series Analysis
A time series analysis describes the changing characteristics of time series in the past and predicts the trends of the future time series; this is an important method for studying the dynamic characteristics, periodic characteristics, and correlation of statistical indicators. Scholars usually use historical data on the incidence rate of vibriosis to predict future incidence rates. However, the incidence rate is related to historical values and also depends on the values of other relevant variables in the same time series, such as temperature, relative humidity, and other climate variables. The covariance and correlation coefficients can be used to describe the degree of correlation between the two time series variables. A correlation analysis determines the interdependence between two or more variables to measure the degree and direction of the correlation between variables and the internal relationship between variables. Pearson's correlation analysis is commonly used and is calculated as follows given a data pair of two time series variables X and Y: Stronger correlations are represented by larger absolute values of the correlation coefficient (correlation coefficients close to 1 or −1), while weaker correlations are represented by correlation coefficient values that are closer to 0.

Spatio-Temporal Scanning Statistics
Spatio-temporal scanning statistics are a common method to detect epidemic patterns of diseases at any time and space in the field of epidemiology [34]; they can be used to detect an increase in the number of local bacterial foodborne diseases and to verify whether the increase is due to random variation. The basic principle of spatio-temporal scanning statistics is to establish a spatio-temporal two-dimensional cylinder activity window based on geographical coordinates, where each cylinder reflects a possible aggregation area [35]. That is, the geographical area and time period of possible disease outbreaks. The bottom of the cylinder represents the geographical area of detection, and the height of the cylinder represents the detection time.
The expected number of cases in each scan window was calculated, and the log likelihood ratio (LLR) was used to evaluate the abnormal degree of the number of cases in the scan window based on the actual number of cases and population [36]. The relative risk (RR) was used to evaluate the risk of bacterial foodborne diseases in the cluster area [37]. The maximum likelihood ratio of all windows corresponds to the most likely cluster. The equation for determining the LLR is shown in Formula (2), and the RR of the cluster area is estimated using Formula (3) under the assumption of a Poisson distribution.
where C is the total number of cases, C A is the actual number of cases in the window, and µ A is the expected number of cases in the window.

Principal Component Analysis to Determine the Effects of Meteorological Factors
Many studies often use original meteorological data for modeling when exploring the dynamic relationship between meteorological factors and dependent variables. However, complex correlations could exist among meteorological factors, and high multicollinearity may increase the uncertainty of the analysis results. A principal component analysis is a multivariate statistical analysis method that reduces the dimension of multiple related variables to several representative comprehensive variables using linear transformations, and it contains information that is not duplicated [38]. The index is obtained after a comprehensive analysis called the main component. The weight determined by this method is based on the internal structural relationship between the indicators obtained by the data analysis; therefore, it is unaffected by subjective factors. The principal components are independent of each other, which helps reduce the interaction of information.

Vector Autoregressive (VAR) Model
The vector autoregressive model is suitable for multivariable time series systems to analyze the dynamic impact of random disturbances on variable systems and explain the effects of various shocks on variables [39]. The model can determine the dynamic relation-ship between climatic variables and the detection rate of bacterial foodborne diseases in a time series; it does not need to distinguish between exogenous and endogenous variables or add constraints to the model. The basic principle is to consider each endogenous variable in the system as a function of the lag value of all other endogenous variables in the system and construct a model. The mathematical expression of the VAR model (p) is where y t is the k-dimensional endogenous variable vector, x t is the d-dimensional exogenous variable vector, ε t is the random disturbance vector, p is the lag order, T is the number of bacterial foodborne disease case samples, and A 1 . . . A p , and B represents the k × k dimensional coefficient matrix that requires estimation. The construction of the VAR model mainly includes the following steps: Unit root test of sample data: It is necessary to test the unit root of the time series before establishing the model to ensure the stationarity of the time series and avoid the pseudo-regression phenomenon. The augmented Dickey-Fuller test (ADF test) was used to evaluate the unit root [40], and the null hypothesis was that there was a unit root. There is no unit root if the sequence is stable following this test. The t-statistic is under three confidence levels (1%, 5%, and 10%), and the original hypothesis is rejected; otherwise, there is a unit root.
The optimal lag order was determined and a VAR model was constructed: A greater order represents a greater degree of freedom for the model. However, a larger order indicated that more parameters required estimation. This partly affected the validity of the model. Therefore, there should be enough lag terms and a sufficient number of degrees of freedom when determining the optimal lag order. This is mainly based on Akaike information criterion (AIC) and Schwarz Criterion (SC), which state that a model with smaller statistics has a better fitting effect. The stationarity test and fitting effect of the model are good if the eigenvalues of each variable calculated by the model fall in the unit circle.
Granger causality test: The causal relationship between variables was examined by testing whether the lagging term of one variable has an impact on the other variables.
Impulse response: The impulse response function is mainly used to analyze the dynamic effects of the random disturbance of each endogenous variable on itself and all other endogenous variables.

Multivariate Time Series Correlation Analysis
A total of 182,473 samples of bacterial foodborne diseases were collected from 2014 to 2018 in the Zhejiang Province. There were 6430 (3.52%) positive cases due to V. parahaemolyticus infection, and 6226 of 6430 positive cases occurred between the months of May and October during these 4 years. Table 1 summarizes the descriptive statistics for the weekly detection rate of V. parahaemolyticus and the ambient meteorological conditions on the detection date. Time sequence diagrams were used to visually show temporal trends in pathogen infection and meteorological variations ( Figure 3). The V. parahaemolyticus infection showed an annual peak during the summer that was consistent with the temperature changes. Moreover, the correlation analysis showed that V. parahaemolyticus was closely correlated with all seven climatic factors, but showed a time lag (Figure 4). Six climatic factors had a positive lag effect on pathogen detection, with different lag periods. The duration of sunlight was over two lag weeks, the temperature was over three lag weeks, and the relative humidity and precipitation were over eight weeks. Meanwhile, wind speed (lag of 11 weeks) had a weak negative correlation (−0.171) with vibriosis.      The PCA model was established based on the seven meteorological indicators X1, X2, …, X7 to reduce the interference of multicollinearity. The cumulative variance contribution rate of the three principal components (F1, F2, and F3) was 90.377%, which meets the requirement of a cumulative variance of 0.85. A further quantitative analysis of the climatic factors showed that the first principal component reflecting temperature-related climate parameters (including mean temperature, daily maximum temperature, and daily minimum temperature) had a variance contribution rate of 46.454%. The second principal component reflects the climate parameters related to water (including average relative humidity, precipitation, and sunlight hours) and had a variance contribution rate of 29.505%. The third principal component was the wind-related climate parameters (including the average wind speed), with a variance contribution rate of 14.418%. The relationships between the principal components and original climate variables are shown in Table  2.  The PCA model was established based on the seven meteorological indicators X 1 , X 2 , . . . , X 7 to reduce the interference of multicollinearity. The cumulative variance contribution rate of the three principal components (F 1 , F 2 , and F 3 ) was 90.377%, which meets the requirement of a cumulative variance of 0.85. A further quantitative analysis of the climatic factors showed that the first principal component reflecting temperature-related climate parameters (including mean temperature, daily maximum temperature, and daily minimum temperature) had a variance contribution rate of 46.454%. The second principal component reflects the climate parameters related to water (including average relative humidity, precipitation, and sunlight hours) and had a variance contribution rate of 29.505%. The third principal component was the wind-related climate parameters (including the average wind speed), with a variance contribution rate of 14.418%. The relationships between the principal components and original climate variables are shown in Table 2. Three significant spatio-temporal aggregation areas were detected using spatio-temporal scanning statistics (Table 3). Vibriosis had an evident clustering tendency in time and space. The period from the 28th to the 37th week of 2016 (June-August 2016) had the highest incidence of vibriosis. The detection rate of V. parahaemolyticus in the eastern coastal areas and northwest Zhejiang Province was significantly higher than that in the southwest area from a spatial perspective ( Figure 5). Three significant spatio-temporal aggregation areas were detected using spatio-temporal scanning statistics (Table 3). Vibriosis had an evident clustering tendency in time and space. The period from the 28th to the 37th week of 2016 (June-August 2016) had the highest incidence of vibriosis. The detection rate of V. parahaemolyticus in the eastern coastal areas and northwest Zhejiang Province was significantly higher than that in the southwest area from a spatial perspective ( Figure 5).

VAR Model Fitting and Prediction
The study area was divided into four regions based on the three clustered areas and one non-clustered area. Four VAR models were constructed to measure the internal relationship between the detection rate of V. parahaemolyticus and the principal components of the three meteorological factors. It is necessary to conduct unit root and cointegration tests before establishing the model to ensure the stationarity of the time series and avoid pseudo-regression. The unit root tests showed that the original sequences of the detection rate of V. parahaemolyticus and the components of meteorological factors in the four

VAR Model Fitting and Prediction
The study area was divided into four regions based on the three clustered areas and one non-clustered area. Four VAR models were constructed to measure the internal relationship between the detection rate of V. parahaemolyticus and the principal components of the three meteorological factors. It is necessary to conduct unit root and cointegration tests before establishing the model to ensure the stationarity of the time series and avoid pseudo-regression. The unit root tests showed that the original sequences of the detection rate of V. parahaemolyticus and the components of meteorological factors in the four regions were a stationary series (Table 4). This was in accordance with the premise of establishing the VAR model. The three meteorological principal components of the four regions were introduced into the four VAR models as explanatory variables. The information of the AIC and SC were combined to determine that the optimal lag periods of C 1 , C 2 , C 3 , and Non_C were 6, 2, 6, and 2 weeks, respectively.
There was a significant one-way Granger causal relationship between the detection rate (C 1 _DR) and the first meteorological principal component (C 1 _F 1 ) or meteorological principal component 3 (C 1 _F 3 ) for C 1 (Table 5). This indicated that changes in climaterelated parameters of temperature and wind speed significantly changed the detection rate of V. parahaemolyticus. Only the first meteorological principal component (C 2 _F 1 ) significantly increased the detection rate for C 2 , C 3 , or Non_C. This indicated that changes in climate-related parameters of temperature could lead to corresponding changes in the detection rate of V. parahaemolyticus in these areas. All four models treated the local detection rate of V. parahaemolyticus as the dependent variable and the principal components of the three meteorological factors as the explanatory variables to form the regression function. The formulae were as follows: In general, all four VAR models performed well, with adjusted R 2 values ranging from 0.573 to 0.727 (Table 6). All unit roots of the four models were in the unit circle. This indicated that the structure of the model was stable.

Impulse Response Analysis of V. parahaemolyticus Detection
The impulse response analysis provided insight into the dynamic reactions between meteorological parameters and the detection rate of V. parahaemolyticus [41]. Figure 6 illustrated that the impulse response function amounts to one standard deviation of the endogenous variable. In the first column, the response of detection rate to the temperaturerelated principal component (F 1 ) shock showed similar performance in the C 1 , C 2 , and C 3 clusters, while the magnitude was slightly smaller in the Non_C region. The responses to temperature-related component shocks were all positive, indicating that temperature significantly increased the detection rate. This promoting effect showed a significant increase in the first 3 weeks, then slowly increased and stabilized, and finally gradually weakened. The positive effect of ambient temperature on the detection rate of V. parahaemolyticus was more significant in the range of C 1 ; the continuous growth time was longer (up to 7 weeks). The detection rate of V. parahaemolyticus slightly decreased in C 3 from 3 to 4 weeks, and the overall trend was like that observed in C 2 . In other words, it gradually stabilized after 3 weeks.
The impulse responses to the moisture-related principal component (F 2 ) shock slightly fluctuated around 0 in C 1 , C 2 , and C 3 in the first 7 weeks, which was consistent with the results of the Granger causality test with lag 2 and 6 ( Table 5); they exhibited a slightly positive effect after 7 weeks. The moisture variable was exogenous in the model of the Non_C region; therefore, it did not participate in the impulse response test.
Similarly, the magnitude of the responses to the wind-related principal component (F 3 ) shock was smaller than F 1 . To a shock of the wind-related principal component, the response of the detection rate of C 1 positively fluctuated until the 6th week; however, the persistence was the smallest among all meteorological principal components.

Discussion
The risk of vibriosis in Zhejiang Province from 2014 to 2018 was significantly aggregated in the spatial and temporal dimensions. The incidence was high from June to August, and the detection rate of V. parahaemolyticus in the eastern coastal areas and northwest was significantly higher than in the southwest. This may be related to the high-temperature and humid environment in the Zhejiang Province during this period. This climate condition is conducive to the growth and reproduction of microorganisms such as V. parahaemolyticus. Thus, food is prone to spoilage. Differences in regional eating habits and economic conditions may be the key factors leading to spatial differences in bacterial foodborne diseases. The three high agglomeration areas of C1, C2, and C3 belong to areas with relatively high per capita disposable income (typical representatives of Ningbo, Hangzhou, and Wenzhou). The residents' ability to consume food outside the home is relatively high. Aquatic products and seafood were the first food species that were suspected to cause V. parahaemolyticus infection [42]. The eastern region (C1 and C2) is close to the port and wharf and is rich in aquatic products. Local residents live near the sea and enjoy eating seafood and aquatic products. This greatly increases the infection rate of V. parahaemolyticus. Therefore, food safety inspection departments in coastal areas need to standardize and strengthen the routine hygiene testing of aquatic products and seafood in the market, especially during the high-incidence summer period. The market supervision department needs to strengthen the hygiene supervision of various restaurants and the training of employees in standard operations. In addition, residents in coastal areas require public

Discussion
The risk of vibriosis in Zhejiang Province from 2014 to 2018 was significantly aggregated in the spatial and temporal dimensions. The incidence was high from June to August, and the detection rate of V. parahaemolyticus in the eastern coastal areas and northwest was significantly higher than in the southwest. This may be related to the high-temperature and humid environment in the Zhejiang Province during this period. This climate condition is conducive to the growth and reproduction of microorganisms such as V. parahaemolyticus. Thus, food is prone to spoilage. Differences in regional eating habits and economic conditions may be the key factors leading to spatial differences in bacterial foodborne diseases. The three high agglomeration areas of C 1 , C 2 , and C 3 belong to areas with relatively high per capita disposable income (typical representatives of Ningbo, Hangzhou, and Wenzhou). The residents' ability to consume food outside the home is relatively high. Aquatic products and seafood were the first food species that were suspected to cause V. parahaemolyticus infection [42]. The eastern region (C 1 and C 2 ) is close to the port and wharf and is rich in aquatic products. Local residents live near the sea and enjoy eating seafood and aquatic products. This greatly increases the infection rate of V. parahaemolyticus. Therefore, food safety inspection departments in coastal areas need to standardize and strengthen the routine hygiene testing of aquatic products and seafood in the market, especially during the high-incidence summer period. The market supervision department needs to strengthen the hygiene supervision of various restaurants and the training of employees in standard operations. In addition, residents in coastal areas require public awareness and education on food hygiene and the prevention of bacterial foodborne diseases.
The detection rate of V. parahaemolyticus in Zhejiang Province may have resulted from the comprehensive effect of climate components based on the PCA and VAR models in this study. Specifically, the occurrence of this bacterial foodborne disease is directly and critically affected by temperature. The influence of moisture variables such as precipitation, relative humidity, and sunlight hours was more complex, and the effect on the model was not prominent. In addition, the comprehensive effect of climate variables had a lagging effect on the infection of bacterial foodborne diseases, and the lagging effect varied between regions. These two findings are consistent with those of Lake et al. in England and Wales [43] and Zhang et al. in Australia [14]. The impulse response function analysis of the four VAR models showed that the temperature meteorological parameters had a positive effect on vibriosis in all four areas of Zhejiang Province. The impact gradually increased and stabilized over time. Previous studies showed a positive correlation between temperature and the occurrence of bacterial foodborne diseases [44][45][46]. In other words, a higher-temperature environment promotes the reproduction of V. parahaemolyticus and increases the probability of food infection, resulting in a positive impact on the infection rate of the disease. Therefore, the detection rate of V. parahaemolyticus significantly increased in the early stage of temperature change. The impulse response results showed that the linkage relationship between the temperature and detection rate varied in different regions. In general, it is difficult for people to promptly respond to sudden changes in temperature, which increases the risk of infection. However, it will raise their awareness of bacterial foodborne diseases and food deterioration and may reduce the occurrence of disease infection once people adapt to the environment, such as in the case of long-term heat waves. The modeling showed that there were regional differences in the response characteristics of bacterial foodborne diseases to temperature variables. The degree of responses of the infection probability to temperature variables in the aggregation areas were significantly higher than that of the non-cluster area, with a fast growth rate and long duration. The impact of temperature on V. parahaemolyticus infection in C 1 was characterized by a long duration (approximately 7 weeks) and strong response intensity (the maximum was approximately 0.011). C 2 was characterized by a large response intensity (the maximum was approximately 0.0125), whereas C 3 was relatively weak (the maximum was approximately 0.008). The response intensity of Non_C to temperature infection in the non-clustered area of bacterial foodborne diseases was significantly lower than the average value (the maximum value was approximately 0.004). These results showed that the sub-regional model avoided the mutual interference of various factors between regions and maximized the comprehensive effect of meteorological factors on the detection rate of V. parahaemolyticus.
The positive response to wind-related principal component shock showed that wind speed can affect the reproduction of V. parahaemolyticus, viability, and time in the environment. Weather conditions with high wind speeds may be conducive to the spread of bacterial pollutants and may indirectly affect the occurrence of Vibriosis; this is consistent with previous findings [47]. The impulse response of detection rate of V. parahaemolyticus to the moisture-related principal components was not obvious in the initial stage. The growth trend of the impulse response curve showed that the moisture meteorological parameters had a lagging effect on vibriosis after being influenced for a period of 7 weeks. Precipitation in the time series diagram showed an obvious lag effect, and the lag response in different regions was different (Figure 7). For example, the detection peak of V. parahaemolyticus in Zhejiang Province corresponded to the detection peak of V. parahaemolyticus after 3-7 weeks. This was consistent with the impulse response of the moisture meteorological parameters in this model. However, it differs from previous studies. For example, there is a significant positive correlation between vibriosis and relative humidity in Korea [21]. Meanwhile, changes in the meteorological conditions (such as humidity) may change the characteristics of survival and transmission patterns of microorganisms such as V. parahaemolyticus, increasing the number of cases of vibriosis [48]. This difference may be due to those studies failing to consider the combined effects of meteorological factors and spatial heterogeneity. These results indicated spatial differences in the effects of meteorological factors on the detection rate of bacterial foodborne diseases. This reflected the effectiveness of the innovative idea of dividing the study area with spatio-temporal scanning statistics to explore the dynamic relationship between detection rates and meteorological factors. parahaemolyticus, increasing the number of cases of vibriosis [48]. This difference may be due to those studies failing to consider the combined effects of meteorological factors and spatial heterogeneity. These results indicated spatial differences in the effects of meteorological factors on the detection rate of bacterial foodborne diseases. This reflected the effectiveness of the innovative idea of dividing the study area with spatio-temporal scanning statistics to explore the dynamic relationship between detection rates and meteorological factors.

Conclusions
This study proposed a framework to evaluate the dynamic relationship between different meteorological factors and the risk of bacterial foodborne diseases based on the spatio-temporal heterogeneity of vibriosis in Zhejiang Province. A VAR model was combined with foodborne diseases to analyze the spatio-temporal risk characteristics and potential climate risk factors of bacterial foodborne diseases.
Vibriosis exhibited significant temporal and spatial aggregation and peaked in the summer. The detection rate of V. parahaemolyticus was significantly higher in the eastern coastal areas and northwest compared with the southwest. In addition, temperature, relative humidity, precipitation, sunlight hours, and wind speed were important meteorological indexes affecting vibriosis that lag by 3, 8, 8, 2, and 11 weeks, respectively. The comprehensive dynamic influence of the meteorological factors showed regional variation. The temperature-related principal component had a stronger promoting effect on the detection rate in spatio-temporal cluster areas, while the influence in the non-clustered area was relatively weak. The zoning modeling framework proposed in this study considered the temporal aggregation effect of vibriosis and the spatial agglomeration characteristics of epidemics from the spatial dimension compared with whole region modeling as conducted in previous studies; this maximized the impact of reduced climate variables on regional foodborne disease infection probability. This study provided guidance and

Conclusions
This study proposed a framework to evaluate the dynamic relationship between different meteorological factors and the risk of bacterial foodborne diseases based on the spatio-temporal heterogeneity of vibriosis in Zhejiang Province. A VAR model was combined with foodborne diseases to analyze the spatio-temporal risk characteristics and potential climate risk factors of bacterial foodborne diseases.
Vibriosis exhibited significant temporal and spatial aggregation and peaked in the summer. The detection rate of V. parahaemolyticus was significantly higher in the eastern coastal areas and northwest compared with the southwest. In addition, temperature, relative humidity, precipitation, sunlight hours, and wind speed were important meteorological indexes affecting vibriosis that lag by 3, 8, 8, 2, and 11 weeks, respectively. The comprehensive dynamic influence of the meteorological factors showed regional variation. The temperature-related principal component had a stronger promoting effect on the detection rate in spatio-temporal cluster areas, while the influence in the non-clustered area was relatively weak. The zoning modeling framework proposed in this study considered the temporal aggregation effect of vibriosis and the spatial agglomeration characteristics of epidemics from the spatial dimension compared with whole region modeling as conducted in previous studies; this maximized the impact of reduced climate variables on regional foodborne disease infection probability. This study provided guidance and geographical support for relevant government departments to prevent and control foodborne diseases in Zhejiang Province. This research framework can be extended to other problems caused by climate and environmental change.
The following limitations and future perspectives were proposed. First, the modeling area in this study was based on the administrative region and was divided according to the results of spatio-temporal scanning statistics. This method ignored differences in social and economic conditions and may be inadequate. For example, the 2017 per capita GDP in Haining City and Taishun County was CNY 103,220 and CNY 36,420, respectively. The gross domestic product of the two places was quite different, but it was divided into non-clustered areas (Non_C) in the spatio-temporal scanning results. An understanding of regional divisions, conducted in a more reasonable manner, will play an important role in future studies. Second, the evaluation framework for the drivers of foodborne disease incidence must be further improved. This study only considered meteorological factors in the modeling; thus, future studies should include indicators such as food consumption structure, food exposure information, age structure, and financial and health expenditure to explore the factors influencing the incidence of bacterial foodborne diseases under spatio-temporal heterogeneity from a more comprehensive point of view.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.