Introduction
Bacteria of the genus Vibrio are found in surface waters in the Americas, Asia, Europe and Australia and are responsible for causing several types of human infections, including cholera and non-cholera vibriosis (NCV) [Reference Armstrong, Hollingsworth, Morris, Evans and Brachman1, Reference Franca2]. While cholera is typically caused by toxigenic Vibrio cholerae serogroups O1 or O139, vibriosis is caused by non-toxigenic V. cholerae serogroups and roughly 12 other Vibrio species including V. vulnificus, V. parahaemolyticus and V. fluvialis [Reference Newton3, Reference Mahon4]. While most individuals with vibriosis develop symptoms such as diarrhoea and abdominal pain, some (particularly those with liver disease or immune compromise) may experience skin and soft tissue infection (typically characterised by bullous skin lesions) and septic shock [Reference Dechet5]. Infection can be transmitted by food or by exposure to contaminated salt-water, and the most severely affected individuals often have immune compromise or liver disease [Reference Dechet5]. Infections display summertime seasonality, which is believed to reflect the thermophilic nature of these bacteria, as well as increased leisure-related salt-water exposure and harvesting of seafood in the summer [Reference Armstrong, Hollingsworth, Morris, Evans and Brachman1].
There has been an increase in NCV reported in recent decades in the United States [Reference Newton3]. In 1988, the US Centers for Disease Control and Prevention (CDC) established the Cholera and Other Vibrio Illness Surveillance (COVIS) system, which initially focused on states (Texas, Alabama and Louisiana) forming part of the Gulf of Mexico coastline and expanded through the 1990s. By 2007, it was established as a national surveillance system at which time NCV became nationally notifiable [Reference Newton3, Reference Wong6]. The system relies on reports of laboratory-confirmed vibriosis cases from state and territorial public health authorities [Reference Newton3, Reference Wong6]. The reason for increased incidence of NCV in the United States is unclear, but might reflect improved case ascertainment, emergence of novel vibrio strains [Reference Martinez-Urtaza7] or increased abundance of the pathogen due to increasing ocean temperatures. Recent, concurrent, increases in NCV in the Baltic region of Europe appear to be linked to warming Baltic Sea-surface temperatures and may be a result of anthropogenic climate change [Reference Semenza8].
Large-scale climatic phenomena may provide an important natural experiment that improves understanding of the potential future impacts of climate change on infectious disease risk. For example, the impact of the El Niño Southern Oscillation (ENSO) on cholera risk has been well-studied [Reference Hashizume9, Reference Pascual10]. The irregular and periodic nature of ENSO, and the similarity of ENSO-linked weather anomalies to those projected to occur with greater frequency with climate change, make ENSO a powerful natural experimental system for understanding the potential impact of climate change on infectious diseases and other health and economic issues [Reference Fisman, Tuite and Brown11]. Oceanic warming and extreme weather events including heat waves and hurricanes have also been associated with range expansion for vibriosis risk [Reference Baker-Austin12] and NCV outbreaks, respectively [Reference Heilpern and Borg13]. To our knowledge, the impact of ENSO and other large-scale climatic phenomena on vibriosis risk in the United States has not been previously examined.
We combined environmental exposure data with national and state-level vibriosis case counts derived from the COVIS system to evaluate the degree to which variation in vibriosis risk in the United States and in U.S. regions might be explained by large-scale climatic phenomena such as ENSO. We also performed exploratory regional analyses to evaluate regional differences in evolving vibriosis risk, and to evaluate the possibility that changing patterns of risk are consistent with vibriosis range expansion.
Methods
Non-cholera vibriosis data
Annual COVIS reports containing vibriosis case counts in the United States from 1997 to 2014 are publicly available, and were obtained from the United States Centre for Disease Control and Prevention [14]. The initial report, which included aggregate data for 1997 and 1998, was not used. Annual reported vibriosis case counts, by state, were extracted manually from published maps. Monthly national case counts were available from published graphs; we used the WebPlotDigitizer application (version 4.1) to accurately estimate case count numbers derived from these figures [Reference Rohatgi15]. To ensure accuracy of extractions, monthly case count estimates were summed and compared with the reported annual national case counts. Incidence was estimated using U.S. national and state-level population denominators obtained from the U.S. Census Bureau [16]. Where necessary, estimates for inter-censal years were estimated via linear interpolation and extrapolation.
Environmental exposure data
As noted above, ENSO is a complex climatic phenomenon manifested by changes in sea-surface temperature and atmospheric pressure over the Pacific Ocean. ENSO has multiple attributes, including sea-surface temperature, wind, air temperature and cloud anomalies. For simplicity, we used the National Oceanic and Atmospheric Administration (NOAA) multivariable ENSO index (MEI) [Reference Wolter17] as our measure of ENSO activity. The MEI is scaled from −3 (consistent with ‘La Niña’-like conditions) to +3, with high values corresponding to ENSO events. Monthly MEI values for January 1999 to December 2014 were obtained from the NOAA [Reference Wolter17].
Although the effects of ENSO on North American weather are widespread, the strongest effects are observed in western states, and to some extent on Gulf states. Weather and ocean temperatures on the east coast of the United States are influenced by the North Atlantic Oscillation (NAO) index [18], which manifests as fairly regular fluctuations in atmospheric pressure between Iceland and the Azores. The resultant pressure gradient influences land and sea-surface temperatures, and precipitation patterns, in North America and Europe. Our measure of NAO strength, which we also obtained from the NOAA, was the NAO Index. Like the MEI, this index is scaled from −3 to +3, with higher values associated with warmer temperatures and increased precipitation in the eastern United States and lower values associated with cooler, dryer conditions [18]. Unlike ENSO, which is defined in part by sea-surface temperature anomalies, the relationship between sea-surface temperature and NAO appears to be more complex. While certain analyses suggest an inverse relationship between NAO and sea-surface temperatures at prolonged (4 year) lags [Reference Xu19], others suggest reversed causation, with NAO linked to preceding sea-surface temperature anomalies in the Gulf Stream [Reference Wang20].
Statistical analysis
National analyses
We evaluated year-on-year and seasonal trends in overall vibriosis incidence in the United States using count-based regression models that incorporated linear, quadratic and cubic multi-annual trends, as well as fast Fourier transforms to account for seasonal oscillation. As both deviance and Pearson's goodness-of-fit statistics suggested overdispersion of count data, negative binomial models were constructed. Census population estimates were used as model offsets with interpolation used to generate populations in inter-censal years.
Initial estimates of the lagged impact of ENSO and NAO on disease risk were generated by incorporating monthly MEI and NAO index values into trend models at lags of 1–12 months. Models were constructed by backward elimination, with covariates retained for P < 0.2. Combined incidence rate ratios (IRRs) and standard errors for ENSO and NAO exposures at multiple lags were generated via linear combination [Reference Hosmer, Lemeshow, Shewhart and Wilks21].
The inclusion of multiple-lagged exposures can result in several challenges, including difficulties with model interpretation, correlation in exposures at multiple lags and non-parsimony of models. We therefore constructed distributed lag non-linear models to evaluate the integrated effects of environmental exposures at multiple lags [Reference Gasparrini, Armstrong and Kenward22]. These models characterise exposures as two-dimensional risk planes, referred to as ‘cross-bases’, and are defined as level of exposure by lag. We created two cross-bases, for MEI and NAO, at values ranging from −3 to +3, across 12 month lags, which were incorporated into generalised linear models that also adjusted for seasonality, and linear and cubic trends, using a log-link function, in order to approximate negative binomial models constructed above. We assumed that ENSO and NAO effects on risk would be linear at any given lag, but modelled lag structure as a cubic polynomial. Both ENSO and NAO effects were evaluated at the uppermost extreme (+3) value, with index values of zero serving as referents.
Regional analyses
COVIS reports annual vibriosis counts by state, with states grouped together according to the nature of their coastlines (Pacific, Gulf, Atlantic or non-coastal). We evaluated relative incidence in each grouping using negative binomial models with state populations used as model offsets. We also evaluated linear trend terms from 1999 to 2014 for each region, and the associations between vibriosis risk and annual mean ENSO and NAO values. The non-coastal region was used as a referent. Regional differences in trend terms, and the effect of average annual ENSO and NAO values, were evaluated using the Cochrane's Q-statistic, with inverse variance weighting [Reference Serghiou and Goodman23].
Linear trends in risk for each state were generated by estimating the year-on-year IRR for vibriosis by state. Negative binomial models for six states (Colorado, Iowa, Nebraska, New Mexico, South Carolina and West Virginia) failed to converge. To evaluate whether trends might exhibit a north-south gradient suggestive of climate change effects, we constructed meta-regressive models using the latitude of state centroids as the explanatory variable [Reference Tuite, Greer and Fisman24]. Our hypothesis was that relative rate of increase in risk due to climate change would be greater at northern than southern latitudes, particularly in coastal states, due to ocean warming.
All data used in this study are pre-collected aggregate counts and are publicly available. Negative binomial models and meta-regression models were constructed using Stata Intercooled version 15 (Stata Corp., College Station, Texas), while distributed lag non-linear models were constructed using the dlnm package for R, version 3.1.5 [25]. Datasets used for analyses can be obtained via Figshare at https://figshare.com/articles/Vibriosis_data_files/6856427.
Results
From 1999 to 2014, a total of 10 800 cases of vibriosis were contained in COVIS reports, with an additional 62 cases noted as part of a V. parahaemolyticus outbreak in Alaska in 2004. Of these 10 862 cases, we were able to assign 10 102 (93.0%) to a month of year in national-level data, while 10 857 (>99.9%) could be assigned to a specific state. Crude national incidence increased threefold during the period under observation, from 0.11 cases per 100 000 population in 1999 to 0.36 cases per 100 000 population in 2014, representing a 7% average annual increase in incidence (IRR 1.071, 95% confidence interval (CI) 1.061–1.081). Significant seasonal oscillation was observed (P < 0.001) with a summertime predominance. As quadratic trend terms were not significantly associated with risk, and cubic trend terms were associated with a worsening of Akaike's information criterion, the final trend model included only a linear term for year, as well as a fast Fourier transform (Table 1).
Temporal trends and impact of environmental exposures on national vibriosis incidence
MEI, multivariable ENSO index; NAO, North Atlantic Oscillation index.
a Per 1-point increase in monthly MEI or NAO.
ENSO was associated with enhanced risk of vibriosis at 3 and 10 month lags while ENSO at a 5 month lag was associated with attenuated risk. The overall effect of combined ENSO coefficients significantly increased vibriosis risk (combined IRR 1.145, 95% CI 1.064–1.231, P < 0.001). NAO effects were retained in the model at five different lags (3, 5 and 10–12 months), and were inconsistent in effect. Overall, the combined effect of NAO was non-significant (IRR 1.064, 95% CI 0.983–1.151, P = 0.12) (Table 1, Fig. 1). Distributed lag non-linear models of ENSO effect were also associated with a significant increase in vibriosis across lags of 1–12 months (integrated IRR 1.940, 95% CI 1.298–2.901). The integrated effect of NAO was again non-significant in distributed lag models (IRR 0.988, 95% CI 0.607–1.607) (Fig. 2).
Reported and model-predicted monthly vibriosis counts, United States 1999–2014. Observed counts are represented by circles; solid curve represents predictions from negative binomial model incorporating fast Fourier transform and year terms. Date is plotted on the X-axis; counts are plotted on the Y-axis.
Association between lagged MEI and vibriosis risk, United States 1999–2014. (a) Risk surface represents modelled association between the MEI (scaled from −3, most La Niña-like, to +3, most El Niño-like, X-axis) and monthly vibriosis risk over 1–12 month lags (Y-axis). Associated RRs are plotted on the Z-axis. (b) Cross-sectional RR associated with the MEI of +3; lagged El Niño-like conditions are associated with downstream integrated RR of vibriosis (RR 1.940, 95% CI 1.298–2.900).
As expected, regional analyses identified highest risks of infection in the Pacific (IRR 13.412, 95% CI 10.607–16.959) and Gulf (IRR 7.752, 95% CI 6.167–9.744) regions, with elevated risk in the Atlantic region (IRR 4.125, 95% CI 3.491–4.873). Overall, yearly average ENSO was associated with a significant increase in vibriosis risk (IRR per 1 unit increase in mean MEI 1.169, 95% CI 1.012–1.350), while yearly average NAO was not associated with increased risk (Table 2).
Temporal trends, regional differences and impact of environmental exposures on annual vibriosis incidence in states
MEI, multivariable ENSO index; NAO, North Atlantic Oscillation index.
a Per 1-point increase in yearly average MEI or NAO.
Significant increases in vibriosis risk were observed in non-Pacific regions during the period under observation, with significant heterogeneity in rate of change across regions (P for heterogeneity <0.001); in the Pacific region, a non-significant increasing trend was observed (P = 0.067). The most marked increases in vibriosis risk were detected in the Atlantic and non-coastal regions (i.e. regions with the lowest incidence at baseline). In regional models, a 1 unit increase in mean MEI was associated with a 60% increase in vibriosis risk in the Pacific region, but no significant effect of ENSO on risk was observed in other regions. No significant change in risk was observed with NAO in any region. Neither ENSO nor NAO effects were associated with significant regional heterogeneity (P for heterogeneity = 0.50 for ENSO; P = 0.52 for NAO) (Table 3 and Fig. 3).
Regional models of annual vibriosis incidence
MEI, multivariable ENSO index; NAO, North Atlantic Oscillation index.
a Per 1-point increase in yearly average MEI or NAO.
State-level incidence, trends and environmental influence on annual vibriosis risk in the U.S. states, 1999–2014. (a) Mean annual vibriosis incidence per 100 000 population, by state; (b) IRRs for year-on-year change in vibriosis incidence from negative binomial models, by state; (c) trends in vibriosis incidence from negative binomial models, by state; (d) IRR for vibriosis risk with a 1-unit change in the MEI, by region. Six grey-shaded states are those for which negative binomial models failed to converge.
State-level IRRs were also markedly heterogeneous (P < 0.001). In meta-regression models, increasing latitude was associated with the relative magnitude of increased relative risk (RR per 10° increase in latitude of state centroid 1.059, 95% CI 1.020–1.099). Adjustment for latitude resulted in a 36% reduction in between-study variance (τ 2) [Reference Abubakar26]. There was no significant association between region and rate of change, after adjustment for latitude. We also found no significant interaction between latitude and whether or not a state was coastal (P for multiplicative interaction term = 0.87) (Fig. 4).
Association between latitude and linear trend in vibriosis incidence, United States 1999–2014. Correlation between state latitude and average yearly IRR for vibriosis. Each circle represents a single U.S. state or District of Columbia, with size inversely proportional to variance in IRR estimates, corresponding to the weight assigned to each state. Circles are colour coded according to COVIS regions. Fitted lines represent the association between state latitude and IRR as predicted using univariable meta-regression (relative change in IRR per 10° increase in latitude of state centroid 1.059, 95% CI 1.020–1.099). Note that six states for which negative binomial models did not converge are excluded.
Discussion
Given that toxigenic V. cholerae is a sentinel global pathogen due to its virulence, capacity for genetic reassortment and epidemic potential, the ecology and dynamics of this pathogen have been extensively studied [Reference Pascual10, Reference Lipp, Huq and Colwell27, Reference Gil28]. However, other Vibrio species of public health importance can also cause severe human disease, and their associated disease burden appears to be increasing. An increase in disease incidence may reflect increased pathogen abundance in salt-water environments, as well as an increase in the prevalence of immune compromised states in the human population, as individuals with immune compromise or iron overload states appear to be at greatest risk for severe vibriosis [Reference Newton3]. Due to the thermophilic nature of Vibrio species and their increased abundance in warmer waters with decreased salinity, there is reason to anticipate that Vibrio infection incidence will increase in coming decades as a result of climate change-associated ocean warming [Reference Baker-Austin12].
Using publicly available Vibrios case count data from the United States COVIS system, we evaluated vibriosis trends and their relation with large-scale oceanic phenomena. We also examined the association between rate of increase and latitude, as we hypothesised that relative increases in risk would be greatest at northernmost latitudes where vibriosis has historically been less common. Our observations suggest that climate change-driven range expansion may be in progress for this disease, and the sensitivity of vibriosis incidence to El Niño-like conditions suggests that climate change is likely to drive further increases in vibriosis incidence in the future. After adjusting for long-term trends and seasonal oscillation, we found vibriosis incidence in the United States to be very sensitive to the ENSO. Using distributed lag non-linear models to examine risk effects at multiple lags, we found that vibriosis incidence in the United States is expected to double in the year after strong ENSO-like conditions (MEI = 3). This finding is consistent with existing work on the ecology of cholera in low- and middle-income countries [Reference Hashizume9, Reference Alajo, Nakavuma and Erume29]. The regional analyses suggested these effects are concentrated in the Pacific region, where states are most strongly teleconnected to ENSO and where baseline vibriosis incidence is the highest.
Strong ENSO effects may provide insight into future changes in disease epidemiology that may occur under climate change scenarios such as increases in temperature, precipitation anomalies and ocean warming effects [Reference Fisman, Tuite and Brown11]. As current anthropogenic climate change is a global phenomenon with no recent precedent, it has been noted that the irregular nature of ENSO allows it to serve as a natural experiment for exploring future climate change effects [Reference Fisman, Tuite and Brown11]. The mechanisms whereby ENSO might increase vibriosis risk are likely to be complex; however, conditions like changes in salinity, ocean warming and increased abundance of nutrients and flow of freshwater into oceans are all associated with ENSO [Reference Hasegawa, Ueki and Ando30, Reference Chavez31], as well as with Vibrio abundance [Reference Constantin de Magny32, Reference Raszl33]. It should be noted that, as with V. cholerae, NCV have complex ecology, and their abundance may be impacted by the abundance of the cyanobacteria and dinoflagellates with which they are associated; these microbial populations are also sensitive to changes in nutrients, water temperature and salinity such as those that accompany ENSO [Reference Chavez31, Reference Morales-Ramirez and Brugnoli-Olivera34]. The complexity of Vibrio ecology, shown by differing effects of ENSO at multiple lags due to its impact on different components of the ecosystem, make the distributed lag approach we have applied particularly attractive.
Since weather and sea-surface warming in the Pacific region are more likely to be influenced by ENSO than the Atlantic region of the United States, we also explored the impact of the NAO on vibriosis incidence. We hypothesised that NAO would be associated with regional vibriosis incidence in the Atlantic and possibly the Gulf regions. However, compared with the large effects observed with ENSO at the national level and in the Pacific region, the effects associated with NAO were inconsistent in national analyses, and absent in regional analyses. Given that NAO is predominantly a pressure phenomenon (rather than a phenomenon explicitly related to ocean temperature like ENSO), the absence of clear and consistent links to vibriosis may be unsurprising.
Range expansion of vibriosis has been reported in Northern Europe and appears to correlate with increasing Baltic-Sea temperatures [Reference Baker-Austin12]. We therefore sought to evaluate the possibility that there was a north-south gradient in rate of change in vibriosis risk in the United States. The meta-regression approach we applied here was identical to the approach applied in our earlier work on Lyme disease [Reference Tuite, Greer and Fisman24]. The gradient we observed would be consistent with range expansion for vibriosis risk as a result of warmer oceans at higher latitudes. However, other mechanisms (such as increased reporting in jurisdictions not previously believed to be risk-areas for vibriosis) could produce similar patterns. This is possible since the COVIS system has observed increased coverage and increased case submission over time, particularly once NCV became nationally notifiable in 2007 [Reference Wong6]. However, given that similar trends in rates of vibriosis have been reported using FoodNet data, which involve active rather than passive surveillance, it is unlikely that these trends reflect a reporting artefact [Reference Newton3].
Furthermore, an ocean warming effect would be consistent with the ENSO effects observed nationally, and in the Pacific region, in our initial analyses. Indeed, the substantial face-validity of the ENSO effects we observe supports the use of this irregular climatic phenomenon as a ‘natural experiment’ that can be used to anticipate impacts of future climatic change on infectious diseases caused by pathogens with environmental reservoirs. Given the unprecedented and global nature of anthropogenic climate change, evaluation of effects in a ‘control’ population or control period is not possible [Reference Hsiang, Meng and Cane35]. While year-on-year increases in disease incidence could be attributable to gradual increase in land and sea-surface temperatures, it is difficult to attribute such changes to climate as opposed to other time-varying effects, such changing characteristics of diagnostic tests or surveillance systems [Reference Kontopantelis36]. In this context, creative use of irregular ecological exposure data may provide insights into climate-driven disease risk that inform policy and prevention efforts.
Our analysis is inevitably subject to limitations. As we only have access to case counts by time period and region, we cannot account for changing reporting behaviour as a driver for observed trends. However, it should be noted that we do account for such temporal trends in our models, and the ENSO effects we observe are robust after accounting for trends, regardless of their underlying mechanism. Our lack of access to individual case data means that we cannot perform nuanced subgroup analyses that consider case comorbidity or exposure history. Extracting case counts from graphical plots produced counts similar to those contained in reports, but is likely to have resulted in some random misclassification of these counts; while this may have resulted in diminished statistical power, such misclassification would not have resulted in spurious associations or systematic bias in effects. Similarly, non-differential misclassification of exposure, as may have occurred in regional analyses if we attributed illness to a particular state whereas exposure occurred in a different state, would have biased our results towards the null, making the (strong) effects we observe here lower bound estimates.
In summary, we used publicly available data on vibriosis incidence in the United States to perform the first (to our knowledge) evaluation of the possible contributions of ongoing environmental change to observed increases in vibriosis incidence in North America. Broadly, we found two lines of evidence suggesting that climatic change and ocean warming may be important drivers of observed increases in risk. First, Vibrio risk was strongly associated with lagged ENSO-like oceanic changes, and second, a north-south gradient in relative increase in risk. Attributing such changes to ocean warming and climatic change is biologically plausible, and consistent with effects observed in other geographic locales with V. cholerae. Our work suggests that vibriosis may be an important sentinel for infectious disease impacts of climate change, and our methods also present a model that can be applied to other infectious diseases. Anticipating future trends in vibriosis as oceans warm may help with projections of likely future burden of disease, and may help guide and prioritise preventive strategies.
Conflict of interest
None of the authors has any conflict of interest, financial or otherwise, associated with this work.
Open access
