The Seasonal Microbial Ecology of Plankton and Plankton-Associated Vibrio parahaemolyticus in the Northeast United States

ABSTRACT Microbial ecology studies have proven to be important resources for improving infectious disease response and outbreak prevention. Vibrio parahaemolyticus is an ongoing source of shellfish-borne food illness in the Northeast United States, and there is keen interest in understanding the environmental conditions that coincide with V. parahaemolyticus disease risk, in order to aid harvest management and prevent further illness. Zooplankton and chitinous phytoplankton are associated with V. parahaemolyticus dynamics elsewhere; however, this relationship is undetermined for the Great Bay estuary (GBE), an important emerging shellfish growing region in the Northeast United States. A comprehensive evaluation of the microbial ecology of V. parahaemolyticus associated with plankton was conducted in the GBE using 3 years of data regarding plankton community, nutrient concentration, water quality, and V. parahaemolyticus concentration in plankton. The concentrations of V. parahaemolyticus associated with plankton were highly seasonal, and the highest concentrations of V. parahaemolyticus cultured from zooplankton occurred approximately 1 month before the highest concentrations of V. parahaemolyticus from phytoplankton. The two V. parahaemolyticus peaks corresponded with different water quality variables and a few highly seasonal plankton taxa. Importantly, V. parahaemolyticus concentrations and plankton community dynamics were poorly associated with nutrient concentrations and chlorophyll a, commonly applied proxy variables for assessing ecological health risks and human health risks from harmful plankton and V. parahaemolyticus elsewhere. Together, these statistical associations (or lack thereof) provide valuable insights to characterize the plankton-V. parahaemolyticus dynamic and inform approaches for understanding the potential contribution of plankton to human health risks from V. parahaemolyticus for the Northeast United States. IMPORTANCE The Vibrio-plankton interaction is a focal relationship in Vibrio disease research; however, little is known about this dynamic in the Northeast United States, where V. parahaemolyticus is an established public health issue. We integrated phototactic plankton separation with seasonality analysis to determine the dynamics of the plankton community, water quality, and V. parahaemolyticus concentrations. Distinct bimodal peaks in the seasonal timing of V. parahaemolyticus abundance from phyto- versus zooplankton and differing associations with water quality variables and plankton taxa indicate that monitoring and forecasting approaches should consider the source of exposure when designing predictive methods for V. parahaemolyticus. Helicotheca tamensis has not been previously reported in the GBE. Its detection during this study provides evidence of the changes occurring in the ecology of regional estuaries and potential mechanisms for changes in V. parahaemolyticus populations. The Vibrio monitoring approaches can be translated to aid other areas facing similar public health challenges.

IMPORTANCE The Vibrio-plankton interaction is a focal relationship in Vibrio disease research; however, little is known about this dynamic in the Northeast United States, where V. parahaemolyticus is an established public health issue. We integrated phototactic plankton separation with seasonality analysis to determine the dynamics of the plankton community, water quality, and V. parahaemolyticus concentrations. Distinct bimodal peaks in the seasonal timing of V. parahaemolyticus abundance from phyto-versus zooplankton and differing associations with water quality variables and plankton taxa indicate that monitoring and forecasting approaches should consider the source of exposure when designing predictive methods for V. parahaemolyticus. Helicotheca tamensis has not been previously reported in the GBE. Its detection during this study provides evidence of the changes occurring in the ecology of regional estuaries and potential mechanisms for changes in V. parahaemolyticus populations. The Vibrio monitoring approaches can be translated to aid other areas facing similar public health challenges. role of plankton in V. parahaemolyticus ecology (6,8,12). We hypothesized that a direct study of the plankton community and associated V. parahaemolyticus dynamics in the GBE would yield a more comprehensive and informative understanding of how plankton may affect V. parahaemolyticus in this region.
Little is currently known about plankton community dynamics in the GBE, and this type of study has not been conducted previously. Therefore, the goals of this work were (i) to identify the taxa in the plankton community, (ii) to determine whether individual taxa or overall plankton concentrations covaried with V. parahaemolyticus concentrations, and (iii) to identify water quality variables that covaried with plankton and associated V. parahaemolyticus in order to improve V. parahaemolyticus modeling and risk forecasting in this region. This study is the first comprehensive, multiyear assessment of plankton and concentrations of V. parahaemolyticus from plankton in the GBE and the Northeast United States.

RESULTS
Sampling began in July 2014, while in 2015 and 2016, sampling began when smallcraft vessels could be safely operated dependent on ice-out and seasonal conditions in March (2016) and May (2015). Sampling ended for the same reason in November (2015) or December (2014 and 2016). Thirty-one total sampling events with complete data for V. parahaemolyticus concentrations, plankton community analysis, water quality, and nutrients were conducted. V. parahaemolyticus concentrations and plankton total abundance varied widely overall (Table 1); however, they were consistent between years (see Table S1 in the supplemental material). Water temperature, salinity, pH, turbidity, and chlorophyll a and total dissolved nitrogen (TDN) concentrations were consistent with data in previous reports (6,8,22).
Vibrio parahaemolyticus and plankton abundance. Plankton were present in every sample throughout the study period. However, V. parahaemolyticus (detection limit is a most probable number [MPN] of .0.03/g) was detected in only 54.8% of phytoplankton samples (17/31) and 45.1% of zooplankton samples (14/31). Phytoplankton concentrations (cells/liter) (Fig. 2, shown in green) were highest during the spring bloom and in later summer months. V. parahaemolyticus concentrations from phytoplankton were also elevated during the summer months, but unlike phytoplankton concentrations, this organism was not detected in spring and late fall months. Zooplankton (Fig. 2, shown in blue) were present in the water column year-round but at low levels during spring and fall. V. parahaemolyticus was not detected in zooplankton samples during cooler months, when zooplankton levels in the water column declined.
Thirty-four individual taxa were identified and classified into the following groups: diatoms, dinoflagellates, zooplankton, and "other" (one Haptophyte and one Chrysophyte) ( Table 2). Diatoms were most abundant (relative abundance [RA] = 96.4%) followed by zooplankton (RA = 2.8%) and dinoflagellates (RA = 1.1%). The diatoms Chaetoceros spp., Helicotheca tamensis, Navicula spp., and Skeletonema spp. were classified as abundant, as they had the highest total and relative abundance, accounting for 89.8% of the total abundance over the entire study period. Fragilariopsis spp. and Navicula spp. were present in all 31 samples. The majority (77.8%) of phyto-and zooplankton taxa sampled during this 3-year period had been previously identified in the GBE (Table 2) (23,24).
Seasonality in the plankton community, Vibrio parahaemolyticus, and water quality variables. Each plankton sample contained an average of 14 unique phytoand zooplankton taxa (minimum = 6; maximum = 21) (Fig. 3). Chaetoceros spp. and Helicotheca tamensis were present at high concentrations and dominated the plankton community during summer months (Shannon diversity index [H] , 2). Samples from the spring and fall months were generally more diverse (Shannon H . 2), with the exception of early spring samples, which were predominantly Skeletonema spp. For example, in 2016, Skeletonema sp. concentrations exceeded 7,700 cells/liter in three samples in the spring months, and these samples were comparatively less diverse than other samples from spring months. Nauplii and Navicula spp. were generally present at low levels throughout the study period, except in July 2015, when both were detected at .1,000 cells/liter.
The plankton community was strongly seasonal, and Skeletonema spp., Chaetoceros spp. and Helicotheca tamensis, and Rhizosolenia spp. were differentiated by indicator species analysis (ISA) into distinct spring, summer, and fall assemblages, respectively (Table 3). These taxa were detected at approximately the same time in each year of the study and in association with other taxa that appeared in spring, summer, and fall months. Only Coscinodiscus spp. and Thalassionema spp. were not detected in multiple years of the study, as they were detected sporadically in 2015 and 2016. This was identified in our multiresponse permutation procedure (MRPP) analysis but did not alter the season-specific outcomes.
The season-specific patterns of the plankton community were also identified by seasonality models 1 and 2 (Table 4) in individual taxa but were less pronounced in the overall abundance of phyto-and zooplankton. Helicotheca tamensis (0.63 variance explained), followed by Rhizosolenia spp. (0.55 variance explained), Chaetoceros spp. (0.39 variance explained), and Skeletonema spp. (0.35 variance explained), was the most strongly seasonal and best fit with harmonic regression (model 2). The estimated peak timing of the summer-specific Helicotheca tamensis and Chaetoceros spp. occurred in August around day 223 6 11 and day 217 6 18, respectively. The seasonality of V. parahaemolyticus concentrations from phytoplankton were also well fit by the harmonic regression in model 2, and peak timing was estimated to occur in August at day 224 6 16. Nauplii and zooplankton-associated V. parahaemolyticus were well fit by model 1 and 2, and peak timing was similar (around day 188 6 34 and 182 6 24, respectively). Taxa that were infrequently detected and at low abundance were not well fit by the seasonal variables in model 1 or 2 (see Table S2).
The variables water temperature, salinity, dissolved oxygen (DO), and PO 4 were best fit by harmonic regression in model 2 (.0.65 deviance explained) and peak timing occurred after day 213. NO 3 1 NO 2 was well fit by models 1 and 2, with peak timing around day 200 6 14. Chlorophyll a, pheophytin, turbidity, and the other measured nutrients had little or no seasonality based on poor fits to both models (Table S1; Fig. S2).
Correlation analysis. A cluster of significant positive correlations between V. parahaemolyticus concentrations in phytoplankton, Helicotheca tamensis, Chaetoceros spp., PO 4 , salinity, and water temperature is present in the upper left corner of Fig. 4. Likewise, clustering and significant correlative relationships between V. parahaemolyticus concentrations in zooplankton, photoperiod, nauplii, and copepods are also present in the upper middle section of Fig. 4. DO, NO 3 1 NO 2 , Rhizosolenia spp., and Stephanopyxis spp. produced a separate significant cluster but were negatively correlated with variables that clustered with V. parahaemolyticus concentrations in phyto-or zooplankton. No significant correlations were identified between chlorophyll a, pheophytin, total suspended solids (TSS), particulate nitrogen (PN), and particulate carbon (PC) and the majority of individual plankton taxa or V. parahaemolyticus concentrations.
Seasonal microbial ecology in the GBE. V. parahaemolyticus concentrations from plankton were bimodal (Fig. 5). V. parahaemolyticus phytoplankton were temporally similar to and correlated with Chaetoceros spp., and Helicotheca tamensis. V. parahaemolyticus from zooplankton was strongly correlated and seasonally similar to the individual components of nauplii and copepods. In early spring and late fall, Chaetoceros spp., Helicotheca tamensis, nauplii, and copepods were detected in the water column, but V. parahaemolyticus was not detected in phyto-and zooplankton. The seasonal water quality will be important for understanding these bimodal patterns and differential detection. Both V. parahaemolyticus are statistically associated with water temperature. However, V. parahaemolyticus concentrations in zooplankton (peak timing at day 188 6 34) were highest prior to the warmest water temperatures, and the highest concentrations of V. parahaemolyticus from phytoplankton (peak timing at day 224 6 16) occurred approximately one week after the peak timing of water temperature. NO 3 1 NO 2 , DO, and photoperiod were most strongly related to nauplii and V. parahaemolyticus concentrations from zooplankton in both correlative relationship and seasonal timing. Salinity and PO 4 were most strongly correlated with V. parahaemolyticus concentrations from phytoplankton and individual phytoplankton taxa. These environmental variables covaried positively with water temperature, unlike DO and NO 3 1 NO 2 , though their maxima occurred approximately a month after the maxima of water temperature and were similar in timing to the individual phytoplankton taxa and V. parahaemolyticus from phytoplankton.

DISCUSSION
V. parahaemolyticus concentrations cultured from plankton in the GBE are highly seasonal, are bimodal, and are most closely associated with individual plankton taxa rather than total plankton abundance. V. parahaemolyticus concentrations from phyto-and zooplankton are differentially associated with water temperature, salinity, pH, DO, NO 3 1 NO 2 , and PO 4 , which suggests niche-specific preferences. Importantly, common water quality variables and indicators of plankton abundance or V. parahaemolyticus associated with plankton, such as TDN and chlorophyll a, were not statistically or temporally associated with the dynamics of either group, which provides important insight into why these variables did not contribute well to previous V. parahaemolyticus modeling efforts (8) in this region. Very few recent data on plankton communities and their dynamics in the GBE were available prior to this study, so the results of this work provide an important update for this area, and this study can be built on for continued long-term, robust assessment of plankton-V. parahaemolyticus dynamics and monitoring in this region.
V. parahaemolyticus concentrations were most closely associated with a subset of plankton taxa and this interaction appears to depend heavily on the prevailing water quality conditions. In general, the presence or concentration of plankton overall, nutrients, or chlorophyll a was not a reliable indicator that V. parahaemolyticus was present at detectable levels. V. parahaemolyticus from plankton was detected only in warm weather months, whereas phytoplankton taxa were present in all samples between March to December. Notably, peak plankton-associated V. parahaemolyticus concentrations (0.018 to 21 cells/liter) did not correspond with peak plankton concentrations (10,000 to 35,000 cells/liter), regardless of the taxa. Similarly, V. parahaemolyticus from zooplankton was intermittent and at low abundance, and often V. parahaemolyticus was not cultured from zooplankton samples, even though zooplankton taxa were present in the water column. The majority of plankton taxa observed throughout this study were diatoms, which is consistent with previous reports on the GBE ecosystem (23,24). Only 7 samples exceeded 10,000 cells/liter, which could be considered a low level of productivity relative to other regions that frequently experience concentrations exceeding 10 6 cells/ liter (25,26). Three taxa were the primary drivers for each of the seven highest-concentration samples: Chaetoceros spp. (n = 5), Navicula spp. (n = 1), and Skeletonema spp. (n = 1). Skeletonema spp. were detected in the cooler spring months in the GBE, and elevated concentrations of Skeletonema spp. were detected when sampling could begin safely in early March or April. The single midsummer "bloom" of Navicula spp. in 2015 was unexpected, as Navicula spp. were generally detected year-round at low levels. Chaetoceros spp. dominated most warm-water samples, and their abundance was statistically associated with Helicotheca tamensis.
The taxa that made up the plankton community in this study have been reported around the world from a broad range of environmental conditions. It is also interesting that H. tamensis was absent in prior analyses of the plankton community of the GBE (23,24). The global distribution of these taxa indicates that they are capable of adaptation to a wide range of environmental conditions (27)(28)(29)(30). Since these taxa display such a wide-ranging geographic and ecological distribution, we suggest due caution and consideration before inferring ecological insights or drawing parallels between growing regions based on taxonomic identification.
The pronounced seasonality of climatic conditions in the Northeast United States presents unique analytic challenges (8). In this temperate region, water temperature is the dominant seasonally driven variable, and so the majority of seasonal environmental variables are highly intercorrelated and related to water temperature (31)(32)(33). For this reason, seasonality models using harmonic regression or photoperiod and peak timing were used to overcome this intercorrelation to consider the individual correlative relationships and temporal dynamics of the individual variables.
Harmonic regression analysis enabled the identification of key seasonal features. For example, the highest concentrations of V. parahaemolyticus from phytoplankton occurred just after the warmest water temperatures in late August and early September (day 224 6 16) and were most correlated with and temporally similar to patterns in water temperature, salinity, PO 4 , and H. tamensis. V. parahaemolyticus concentrations from zooplankton, however, were highest earlier in the year around mid-June to early July (day 188 6 34), preceding the warmest water conditions and at approximately the time when nauplii were most concentrated in the water column. Nauplii and V. parahaemolyticus from zooplankton were also positively correlated and temporally associated with photoperiod and inversely with DO, pH, and NO 3 1 NO 2 , i.e., a different set of water quality parameters than those associated with phytoplankton. Work from other regions has also observed that V. parahaemolyticus from different matrices can have distinctly associated water quality conditions (2,33,34). However, this is the first study to identify the temporal and seasonally sequential dynamics of V. parahaemolyticus concentrations between matrices, which can serve as an important focus for future work.
A relationship between photoperiod and plankton community dynamics with V. parahaemolyticus concentrations was also observed by Gilbert et al. (35) and Nilsson et al. (34) in other coastal areas. They suggest that the major driver of microbial community assemblage is likely photoperiod-driven water temperature, and the fine-scale variation and turnover could be attributed to more subtle microcosm-level dynamics like shifting nutrient availability or water quality conditions. In our study, V. parahaemolyticus concentrations from zooplankton and total zooplankton, as well as NO 3 1 NO 2 , were strongly associated with photoperiod, whereas V. parahaemolyticus from phytoplankton, H. tamensis, water temperature, salinity, and PO 4 and peak measurements estimated with harmonic regression occurred approximately 1 month after the longest day of the year, corresponding with earlier suggestions of unique ecological drivers for different components of the microbial community.
Nutrient availability is likely one of the most important limiting factors of plankton and V. parahaemolyticus dynamics (17,34,36). However, the majority of nutrient variables were not temporally or significantly correlated with V. parahaemolyticus in plankton or to individual plankton taxa in this study. Likewise, chlorophyll a is one of the most frequently used variables to study plankton and V. parahaemolyticus dynamics, though the strength of this association is also variable between studies (2,32,37).
In this study, the concentration of chlorophyll a was generally low (6.3 6 4.5 mg/liter). Chlorophyll a was not correlated with either plankton dynamics or taxon blooms, nor was it associated temporally or statistically with V. parahaemolyticus in zoo-or phytoplankton. There is work that has shown that V. parahaemolyticus can also support a free-living lifestyle in the water column by subsisting off nutrient-rich floccules or polysaccharide exudate (35,38). This independent lifestyle could provide insight into why V. parahaemolyticus dynamics are not related to chlorophyll a or observed plankton blooms in the GBE. Alternate variables to chlorophyll a should be considered to monitor plankton dynamics in the GBE, at least at the size fraction of .53 mm.
The dynamics between nutrients, plankton and V. parahaemolyticus are complex, and a lack of significance between variables here could relate to some dimensions of these dynamics that were not accounted for in this study, like the timing of nutrient loading events and both plankton blooms and their decline (39)(40)(41)(42)(43). Future work with more frequent sampling could provide an improved resolution of the contribution of nutrients to the microbial community dynamics in the GBE, but they do not appear to be directly related enough to be useful indicators of plankton-associated V. parahaemolyticus dynamics in the GBE at this temporal resolution.
Conclusions. We developed an in-depth ecological study to describe the plankton-V. parahaemolyticus dynamic in the GBE and assess if the plankton community could contribute to observed in changes V. parahaemolyticus, as it has been suggested for other regions around the world (3). The outcome of this work suggests that the plankton-V. parahaemolyticus dynamic in the GBE is highly seasonal, differs between phyto-and zooplankton, and is dependent on prevailing water quality conditions. Importantly, commonly used variables such as chlorophyll a, most nutrients and overall plankton abundance may not be informative leading indicators of plankton or V. parahaemolyticus dynamics and disease risk in this region. Though it is unclear, at this stage, how this relationship may function long-term, this preliminary study provides a basis for characterizing this complex and important ecological relationship in this region that can be integrated into long-term, ongoing surveillance to more effectively study, monitor, and manage V. parahaemolyticus disease risk for the Northeast United States.

MATERIALS AND METHODS
Study site and environmental sampling. The study area was the Great Bay estuary of New Hampshire, near Nannie Island (NI), which has an important oyster (Crassostrea virginica) bed and is a long-term monitoring location (8,12) (Fig. 1).
Continuous (15-min intervals [Q15]) water temperature, salinity, dissolved oxygen (DO), pH, and turbidity data were obtained from the Stormwater Management Program (SWMP) from 2014 to 2016 for times simultaneous with and preceding sampling events at a SWMP datasonde site (station GRBGB) in close proximity (43.07220, 70.86940) to the NI study site. Grab water samples were collected on each sample date and processed for water quality analyses, including nonpurgeable organic carbon (NPOC), total dissolved nitrogen (TDN), nitrate and nitrite (NO 3 1 NO 2 ,) ammonium (NH 4 ), orthophosphate (PO 4 ), dissolved organic nitrogen (DON), total suspended solids (TSS), particulate carbon (PC), particulate nitrogen (PN), chlorophyll a (CHL) and pheophytin (PHEO) measurements. Data for these parameters were also obtained from the SWMP database (https://cdmo.baruch.sc.edu/dges/). Separate water temperature, salinity, pH, and dissolved oxygen (DO) measurements using YSI 6600 and EXO multiprobe sondes (Yellow Springs Instruments, Yellow Springs, OH), were collected concurrently with sampling at low tide from the NI study site.
Plankton collection and sample processing. Two plankton samples were collected with 53-mm mesh netting. In year 1, a student net (Aquatic Instruments, FL, USA) was used during 10 weighted tows to collect ;140 liters of water. In year 2, a Niskin sampler with a student net was used to collect 160 liters. In year 3, a 30-liter Schindler-Patalas trap (Wildco, FL, USA) with a student net was used to collect 180 liters. All sample counts were transformed to cells per liter to standardize comparisons between years. One sample was phototactically separated using published methods and the plankton separation equipment described previously (44). The separated phyto-and zooplankton fractions were filtered in 53-mm Nitex bolting cloth (Wildco, FL, USA) filter cones. The separated and filtered fractions were air dried and weighed. The second plankton sample was filtered, resuspended to a volume of 50 ml in ,53-mm filtrate water from NI, and preserved with 1% sucrose formalin (45).
Plankton identification and enumeration. The second plankton sample was analyzed with a phase-contrast microscope at Â400 magnification (Olympus, USA) to count 10 nonconsecutive columns (100 quadrants of the Sedgwick rafter grid). Plankton enumeration and concentration determinations were consistent with standard methods for plankton analysis (45) using a phase-contrast microscope and a 1ml grafted Sedgwick rafter (Wildco, FL, USA). Phytoplankton were identified to the genus or species level and zooplankton to higher taxa or functional groups. Plankton identification was confirmed with resources from the work of Dolan and Cooper and Baker et al. (46,47). Sample dilution was conducted, when necessary, for samples that were too abundant for accurate counting using deionized (DI) water.
V. parahaemolyticus concentration enumeration. Phyto-and zooplankton samples were analyzed for associated V. parahaemolyticus concentrations according to previously published methods with a 3-tube MPN alkaline peptone water enrichment and Vibrio CHROMagar (48). Probable V. parahaemolyticus isolates were confirmed by PCR detection of the tlh gene (6,8,48). V. parahaemolyticus concentrations associated with phyto-and zooplankton were transformed from MPN per gram to MPN per liter. Samples without detected V. parahaemolyticus were determined to be below the limit of detection and were assigned a standard value of ,0.03.
Statistical analysis. All statistical analyses were performed in the R Statistical Program and Environment, version 3.5.3 (49), and the vegan Community Ecology Package, version 2.5-2 (50). MPN values for V. parahaemolyticus concentrations were log transformed and plankton counts were log11 transformed to approximate normality and reduce skewness. Significance for all analyses was determined by a P value of ,0.05.
Plankton abundance analysis. Plankton community total abundance (TA), relative abundance (RA), species richness, evenness, and Shannon H diversity were calculated for the study. The taxa observed in collected samples were classified as abundant (.4% of sample total abundance), common (#4% and $0.1%), or rare (,0.1%) in the GBE plankton community. Rare taxa (abundance , 0.1%) were not included in univariate, multivariate, or seasonality analysis.
Seasonality. Seasonal community assemblage was assessed between calendar seasons (spring, March to 21 June; summer, 22 June through 21 September; and fall, 22 September to December) and years (2014, 2015, and 2016) by permutational multivariate analysis of variance (PERMANOVA), multiresponse permutation procedure (MRPP), and indicator species analysis (ISA). The seasonality and general trend for all variables were also explored with models 1 and 2: to fit periodicity and trends of the seasonal oscillations (6). In both models, Y t is the daily time series for the outcome of interest, b 0 is the intercept, t is the daily time series, and b 1 indicates a general trend in the variable of interest. Model 1 contains the photoperiod variable b p photoperiod, and model 2 uses harmonic regression terms for the calendar day in the study where b s and b c are the coefficients of the harmonic terms and v is the term representing the annual cycle (365.25 days, v = 1/365.25). The peak timing of the periodic oscillations identified by model 2 was determined by calculating the phase shift: Whenb S andb c were positive, k was 0. Whenb S was ,0 andb c was .0, k was equal to 2p , and whenb S andb c were negative, or whenb S was .0 andb c was ,0, then k was equal to p . The phase