Impact of hypoxia on gene expression patterns by the human pathogen, Vibrio vulnificus, and bacterial community composition in a North Carolina estuary.

Estuarine environments are continuously being shaped by both natural and anthropogenic sources which directly/indirectly influence the organisms that inhabit these important niches on both individual and community levels. Human infections caused by pathogenic Vibrio species are continuing to rise, and factors associated with global climate change have been suggested to be impacting their abundance and geographical range. Along with temperature, hypoxia has also increased dramatically in the last 40 years, which has led to persistent dead zones worldwide in areas where these infections are increasing. Thus, utilizing membrane diffusion chambers, we investigated the impact of in situ hypoxia on the gene expression of one such bacterium, Vibrio vulnificus, which is an inhabitant of these vulnerable areas worldwide. By coupling these data with multiple abiotic factors, we were able to demonstrate that genes involved in numerous functions, including those involved in virulence, environmental persistence, and stressosome production, were negatively correlated with dissolved oxygen. Furthermore, comparing 16S ribosomal RNA, we found similar overall community compositions during both hypoxia and normoxia. However, unweighted beta diversity analyses revealed that although certain classes of bacteria dominate in both low- and high-oxygen environments, there is the potential for quantitative shifts in lower abundant species, which may be important for effective risk assessment in areas that are becoming increasingly more hypoxic. This study emphasizes the importance of investigating hypoxia as a trigger for gene expression changes by marine Vibrio species and highlights the need for more in depth community analyses during estuarine hypoxia. ©2017. The Authors.

Introduction

Aquatic hypoxia, defined by <30% dissolved oxygen (DO) saturation (<2 mg O2 L−1), in marine systems has been increasing in frequency, severity, and duration by both natural and anthropogenic sources, which is being shaped by global climate change [Conley et al., 2011; Halpern et al., 2008; Rabalais et al., 2002, 2007]. Primary sources of hypoxia are weather patterns such as El Niño and La Niña, coastal upwelling, and anthropogenic eutrophication [Altieri and Gedan, 2015; Diaz and Rosenberg, 2008]. Since the 1970s cases of reported hypoxic zones have increased dramatically with some of the more recent cases being documented in the Baltic Sea, an area undergoing rapid sea surface temperature warming [MacKenzie and Schiedek, 2007; Meier et al., 2016]. Recently, infections caused by such pathogenic Vibrio species as V. vulnificus, V. parahaemolyticus, and V. cholerae are being reported in areas where such infections were never seen before, as far as 65°N, during unseasonably warm heat wave phenomena [Baker‐Austin et al., 2013; Bier et al., 2015; Levy, 2015]. Although the aforementioned increases were not specifically investigated in the context of hypoxia, it has been shown that hypoxia is positively correlated with temperature and is at its peak during summer months worldwide [Rabalais et al., 2009], corresponding to high levels of Vibrio spp. in similar environments [Froelich et al., 2015; Givens et al., 2014; Robles et al., 2013]. Although numerous cases of Vibrio infections have been reported in coastal areas experiencing frequent hypoxic/anoxic events, or from seafood originating from these areas, its impact on virulence potential has been largely overlooked.

A natural inhabitant of estuarine waters worldwide, V. vulnificus, is a potentially fatal human pathogen. This bacterium is responsible for causing septicemia and/or wound infections carrying 50% and 25% mortality rates, respectively [Jones and Oliver, 2009]. V. vulnificus is commonly isolated from a variety of sources, including shellfish, fish, water, and sediment, in estuaries worldwide [Oliver, 2006]. A remarkable feature of this bacterium is its genotypic variability, having not only three distinct biotypes but multiple genotypes [Amaro and Biosca, 1996; Danin‐Poleg et al., 2013; Warner and Oliver, 2008]. Within biotype one strains of V. vulnificus there are two genotypes, clinical (C) and environmental (E), which correlate to both isolation source and specific polymorphisms in the virulence‐correlated gene [Rosche et al., 2005, 2010]. C‐genotype strains are those most frequently causing primary septicemia after the consumption of raw or undercooked seafood, predominantly oysters, whereas E‐genotypes predominate in estuarine environments but appear less virulent [Oliver, 2013; Rosche et al., 2010]. However, there is a subset of E‐genotype strains capable of causing severe wound infections, which generally result from entry into a previously acquired wound or during an injury in estuarine waters harboring V. vulnificus [Oliver, 2005; Ruppert et al., 2004]. Interestingly, the Gulf Coast is home to the largest hypoxic “dead” zone in the U.S. [Diaz and Rosenberg, 2008; Vugia et al., 2013], and 21% of total Vibrio cases in the Gulf Coast states were confirmed to be caused by V. vulnificus. In contrast, non‐Gulf Coast states had only 7% of such cases, as reported by the most recent “Cholera and Other Vibrio Illness Surveillance” Annual Report in 2014 (https://www.cdc.gov/vibrio/surveillance.html).

Oxygen deprivation is experienced by V. vulnificus not only in estuarine waters but also in mollusks during low tide, during postharvesting treatments, and in the gastrointestinal tract of human hosts [Anderson, 2013; Breitburg et al., 2015; He et al., 1999]. It is known that oxygen availability can influence the behavior of bacteria, from rapid changes in metabolic function to the production of virulence factors. Invasion from the GI tract, production of lipopolysaccharide and flagella, and upregulation of key colonization factors are known to be facilitated by exposure to anoxia in vivo by Salmonella, Escherichia coli, and V. cholerae, respectively [Jones et al., 2007; Lamas et al., 2016; Liu et al., 2011; Xu et al., 2003]. We have previously shown that anoxia leads to a translucent colony phenotype and decreased biofilm formation by V. vulnificus, due to the downregulation of genes involved in capsule production and type IV pilus formation in vitro [Phippen and Oliver, 2015a, 2015b]. In situ gene expression by V. vulnificus has been investigated thoroughly in response to both temperature and salinity, and during entry into the viable but nonculturable state [M. K. Jones et al., 2008; Smith and Oliver, 2006a, 2006b]. Such studies revealed that both temperature and salinity impacted the expression of various genes involved in stress response and putative virulence factor production; however, additional biotic and abiotic factors during the sampling events were not considered [M. K. Jones et al., 2008; Smith and Oliver, 2006a, 2006b]. The frequency at which V. vulnificus experiences hypoxia, combined with the lack of understanding which factors are required for virulence, makes it an ideal model for studying both in situ and in vitro hypoxia.

Along with understanding the role of oxygen on gene expression patterns by individual species, it is also important to understand how bacterial communities are impacted. Environmental conditions, such as temperature, salinity, pH, nutrient availability, and oxygen saturation, have been documented to have large impacts on prokaryotic community structure in marine systems. Oxygen availability is a key factor in shaping both fresh and saltwater bacterial community composition, especially when combined with spatial sampling [Bouzat et al., 2013; Spietz et al., 2015]. However, few studies have described these communities as they experience daily estuarine hypoxia, but rather typically focus only on oxygen minimum zones [Stevens and Ulloa, 2008]. Studies examining temporal open‐ocean prokaryotic communities have shown that among other factors, oxygen shows predictability in their structure, suggesting important ecological roles for these species [Fuhrman et al., 2006]. Additionally, Spietz et al. [2015] found that bacterial communities can be impacted at concentrations of DO higher than the definition for hypoxia, which proceeds visible effects on macroorganisms and suggests that these changes occur even in areas experiencing oxygen stress but not hypoxia. This heightened sensitivity to changing oxygen concentrations could potentially explain why certain bacterial species dominate in one condition over another. Furthermore, as global climate change is linked to increasing hypoxic events [Altieri and Gedan, 2015; Bakun et al., 2015; Doney et al., 2012; Justic et al., 1996, 1997], which may impact the distribution of pathogenic species in the environment, the need to understand these communities as a function of environmental change is paramount [K. E. Jones et al., 2008; Vezzulli et al., 2013, 2016].

Here we explore the hypothesis that estuarine hypoxia results in distinct gene expression patterns by V. vulnificus and drives bacterial community composition at a moderately eutrophic site (as designated by NOAA, 1996) in the Bogue Sound of North Carolina. We addressed this hypothesis by (1) monitoring multiple abiotic factors over a 3 day period to determine what parameters were most variable, (2) identifying expression levels of multiple gene families during high‐ and low‐oxygen sampling events, and (3) determining bacterial community compositions during both hypoxia and normoxia.

Methods

Strains and Culture Conditions

A clinical isolate of V. vulnificus, CMCP6, was utilized for this study and was stored at −80°C in Bacto Luria‐Bertani broth (LB) (BD Difco, Franklin Lakes, NJ, USA) with 20% glycerol. This strain was regularly grown overnight (ON) in Bacto heart infusion broth (BD Difco, Franklin Lakes, NJ, USA), at 30°C with aeration in a rotary incubator. For use in chamber deployment, ON cultures were grown and washed three times in sterile 20‰ artificial seawater (ASW: Instant Ocean; Aquarium Systems, Mentor, OH, USA) and resuspended in 250 mL 20‰ ASW.

Site Description and Monitoring of Environmental Parameters

All experiments were conducted in mid‐August at Hoop Pole Creek, which is part of the Hoop Pole Reserve on the eastern North Carolina coast, USA. This site is located near a storm drain input and is frequently exposed to increased nutrients due to rain water runoff. To monitor abiotic parameters, a Hach Hydrolab minisonde 4a data logger (OTT Hydromet, Germany) was utilized for the entirety of the exposure period. Temperature, salinity, pH, and DO were recorded every 30 min for 5 days during the sampling time. The data logger was deployed at morning low tide and was secured to permanent poles positioned in the middle of the Hoop Pole Creek site. Additionally, the data logger was checked twice daily to remove any debris or mud which might impact measurements. Measurements at each sampling event are shown in Table 1.

Table 1

Environmental Parameters Measured at Each Sampling Eventa

Parameter	8/18/15	8/19/15	8/20/15
High DO
DO (%)	83.4	86.4	96.3
Temperature (°C)	30.05	31.69	30.21
Salinity (ppt)	31.69	30.21	29.79
pH	8.13	8.49	8.26
Low DO
DO (%)	21	18.7	22.1
Temperature (°C)	26.19	27.75	28.39
Salinity (ppt)	26.19	29.77	27.54
pH	8.13	7.98	7.59

Measurements recorded with Hach Hydrolab minisonde 4a data loggers.

Membrane Diffusion Chamber Preparation and Deployment

V. vulnificus was grown and prepared as described above and deployed into large volume membrane diffusion chambers. These were modified from the chambers originally designed by McFeters and Stuart [1972]. These modified chambers consist of two 76 mm, 0.2 µm hydrophilic polycarbonate filters (Midland Scientific Inc., Omaha, NE, USA) positioned between two doughnut shaped sections of Plexiglas with a 2 inch piece of polyvinyl chloride (PVC) piping between them (Figure S1 in the supporting information) [McFeters and Stuart, 1972]. Before assembly, PVC pipes were sterilized by a 10% bleach bath for 20 min followed by rinsing in a sterile deionized water bath. Membranes were autoclaved, attached to the PVC pipe by aquarium grade silicon, and each side was checked for leaks by submerging the chamber in water and flowing air through the tube to make sure that no bubbles were produced. Bacterial culture (250 mL) at a final cell concentration of approximately 106 CFU/mL was injected into each chamber using sterile syringes. Once filled, chambers were wrapped in black plastic hex fencing, placed in coolers, and covered with 20‰ ASW for transport to the field site. Chambers were attached securely to the pole harboring the data logger during late afternoon low tide.

In Situ RNA Harvesting

Approximately every 12 h (at each low tide) for three consecutive days duplicate samples were removed from each of three chambers using sterile syringes prefilled with RNAprotect (Qiagen, Valencia, CA) at a 2:1 ratio of RNAprotect to cell culture, following the manufacturer's protocol. A previously optimized RNA extraction assay was then performed as described by Williams et al. [2014a]. RNA quality and quantity was assessed using a NanoDrop spectrophotometer (Thermo Scientific, Wilmington, DE), and samples having a 260/280 nm ratio of ≥1.7 and a concentration of ≥75 ng/μL were stored at −80°C. Using endpoint PCR, determination of DNA contamination was performed as previously described [Phippen and Oliver, 2015b; Williams et al., 2014a] utilizing the species‐specific gene target, vvhA.

Primer Design

Sequenced C‐ and E‐genotype strains of V. vulnificus were used as reference strains for qRT‐PCR primer design, employing NCBI Primer‐BLAST software. Assessment of primer quality and specificity was performed as previously described employing in silico PCR and the IDT OligoAnalyzer 3.1 software [Phippen and Oliver, 2015b; Williams et al., 2014a]. Primer pair efficiencies were predicted by in silico PCR estimation tools, and primers with efficiencies of ≥1.5 were selected and purchased from Sigma‐Aldrich. Primers optimized for this study are listed in Table S3.

Relative qRT‐PCR

A previously described qRT‐PCR protocol was adapted for all in situ expression studies [Williams et al., 2014a]. Gene expression of V. vulnificus strain CMCP6 was examined, in triplicate, for each sampling time using PerfeCTa SYBR green FastMix, Low ROX (Quanta Biosciences, Beverly, MA). Normalization of target genes was performed using glyceraldehyde‐3‐phosphate dehydrogenase, an endogenous control, to correct for sampling error. The Pfaffl equation was used to calculate fold change in gene expression, accounting for differences in the primer pairs, utilizing a previously described PCR efficiency analysis [Pfaffl, 2001]. Gene expression results were analyzed comparing each sample time to the previous sample time, where t = 0 was immediately upon chamber deployment; for example, x = 12 represents fold change of 12 h to 0 h and x = 24 represents fold change of 24 h relative to 12 h. Nonparametric Mann‐Whitney rank‐sum tests were used to determine significance between target transcripts with adjusted p values calculated using the Bonferonni method. All gene expression data were analyzed using GraphPad Prism (version 5.0; GraphPad Software Inc.). Pearson's product‐moment correlational analysis was performed on gene expression and environmental parameters in SPSS Statistics software (version 24, IBM Corp., Armonk, NY, USA) and output visualized using GraphPad Prism.

Water Sampling and DNA Extraction

Water samples were taken at each RNA harvesting event for determination of bacterial community composition. Briefly, 1 L of seawater was removed, stored in sterile 1 L bottles, and placed in a cooler to be transferred to the lab for processing. Water was sampled approximately 3 m from chamber location and sampled at midwater depth. Triplicate water 250 mL samples were passed through a 0.2 µm pore size polyethersulfone filter using vacuum filtration, and the filters kept at −80°C until extraction. DNA from filters was extracted using the FastDNA SPIN kit for soil (MP Biomedicals, Santa Ana, CA, USA) following the manufacturers protocol.

16S rRNA Library Preparation and Analysis

Bacterial sequences were amplified using primers for the V3 and V4 region (Forward: 5′ CTACGGGNGGCWGCAG, Reverse: 5′ GACTACHVGGGTATCTAATCC) as previously described [Klindworth et al., 2012] and triplicates were combined into one sample. Cleanup was performed using the Axygen AxyPrep Kit (Thermo Scientific, Wilmington, DE, USA), followed by index PCR using the Nextera XT Index kit (Illumina Inc., San Diego, CA, USA), and a second PCR cleanup. Fluorometric quantitation was performed after each cleanup step by Qubit (Life Technologies, Carlsbad, CA, USA), and libraries were validated by identification of a ~630 bp product using a Bioanalyzer DNA 1000 chip (Agilent Technologies, Santa Clara, CA, USA). Illumina MiSeq sequencing was performed at the David H. Murdock Research Institute, Kannapolis, NC, USA. Sequences were assembled and edited using PEAR [Zhang et al., 2014] and Trimmomatic 2 [Bolger et al., 2014] and screened for singletons. Using the Qiime pipeline [Caporaso et al., 2010b], chimeric sequences were removed using both ChimeraSlayer [Haas et al., 2011] and USEARCH 6.1 [Edgar et al., 2011] and not used for subsequent analysis of samples. Operational taxonomic units (OTUs) were assigned using MOTHUR [Schloss et al., 2009], aligned using PyNAST [Caporaso et al., 2010a], and assigned to taxonomies using Ribosomal Database Project classifier [Wang et al., 2007] and the Greengenes database [McDonald et al., 2012]. Beta diversity was calculated using both weighted and unweighted UniFrac [Lozupone and Knight, 2005] and visualized using Emperor [Vázquez‐Baeza et al., 2013]. Sequences were made available through the NCBI BioSample database (accession numbers: SAMN05771570, SAMN05771571, SAMN05771572, SAMN05771573, SAMN05771574, and SAMN05771575).

Results and Discussion

Environmental Profiles

In order to determine which environmental parameters might be driving the gene expression patterns of V. vulnificus, we measured DO, temperature, salinity, and pH every 30 min during in situ chamber deployment. Measurements from the datasondes revealed distinct temporal abiotic profiles (Figure 1a). Peaks were designated as “high DO” and valleys as “low DO,” corresponding to trends in DO throughout our sampling periods. As expected, temperature was negatively correlated with DO (Pearson r = −0.917) and ranged from 26.19°C to 31.69°C, whereas pH (7.6–8.5) and salinity (27.9‰–31.7‰) did not correlate with DO (r = 0.684 and r = 0.533). To determine which of these parameters varied the most, we normalized change from one sampling time to the next while removing the unit of measurement, which was key in that each parameter had a different unit and thus could not be statistically compared. Figure 1b shows a significant (one‐way analysis of variance (ANOVA), p < 0.0001) difference in changing oxygen concentrations when compared to change in temperature, salinity, and pH. However, the latter were not statistically significant from one another, further illustrating the highly variable nature of oxygen availability at this site. Average DO saturation during high‐DO sampling times was 88.7%, while low‐DO sampling events fell to an average of 20.9% (Table S1). Interestingly, apart from DO saturation, all parameters were more variable during the low‐DO sampling times but overall remained relatively constant (Table S1).

Figure 1

(a) Temporal environmental parameters as measured every 30 min for three days and (b) normalized change in DO, temperature, salinity, and pH over time. Normalized change was calculated by calculating change from one sampling event to the next followed by subtracting the mean from each measured change and dividing by its SD; values used fell within ±1.5 h around high‐ and low‐DO sampling times. One‐way ANOVA revealed a significant (p < 0.0001) difference in DO normalized change when compared to all other parameters.

Hoop Poll Creek has been described as a highly dynamic site experiencing large daily fluctuations in numerous abiotic factors, especially oxygen [Khan et al., 2012]. This is consistent with what we found for DO; however, the other factors measured did not vary significantly (Figure 1b). Furthermore, this site is potentially eutrophic as a storm drain located at the head of the estuary sporadically adds excess nutrients to the surrounding water, which could exacerbate hypoxic intervals. The combination of the excess nutrients and warm summer temperatures may have resulted in increased photosynthetic activity during the day, resulting in increased DO, and rapid utilization and depletion of this DO overnight [Malone, 1991]. There was no precipitation during our 3 day sampling, which may explain why salinity remained high and relatively consistent. Additionally, there was little variability in daily high temperatures, which was 29.2°C ± 1.2°C. Our results provide further evidence that this site may be used for analyzing the effect of cyclical hypoxia and eutrophication on marine microorganisms.

In Situ Gene Expression Profiles

To characterize the response of V. vulnificus to these cyclical periods of hypoxia and normoxia, duplicate samples were removed from each of the three chambers over 3 days for relative qRT‐PCR analysis during both low‐and high‐DO intervals. Exact environmental parameters measured at each sampling interval are summarized in Table S2. We investigated a variety of genes involved in multiple categories, including capsular polysaccharide (CPS)/exopolysaccharide (EPS) production, motility and attachment, global gene regulation, putative virulence factor production, stressosome formation, and general metabolism (Figure 2). Many of these genes have been previously shown to be differentially regulated by V. vulnificus in response to environmental and host stressors and/or those that have been implicated in virulence [Bisharat et al., 2013; Phippen and Oliver, 2015a, 2015b; Williams et al., 2014b].

Figure 2

Relative temporal expression of representative genes involved in (a) CPS/EPS production, (b) attachment and motility, (c) global gene regulation, (d) putative virulence factor production, (e) stressosome formation, and (f) metabolic function. Temporal fold change represents expression at corresponding x axis time relative to previous sample. Error bars represent SD of three biological replicates (n = 3), comprising two technical replicates averaged for strain CMCP6. Grey shaded areas indicate hypoxic sampling events. Symbols (#/$/) represent values that were not statistically different from the previous time point (one‐way ANOVA, p > 0.05). Arrows indicate that expression was outside the limit of detection at T = 0 in which case the fold changes were set to an arbitrary 100‐fold increase.

Polysaccharide Production, Attachment, and Motility

V. vulnificus produces multiple types of polysaccharides, notably CPS and EPS, which have implications in evading host phagocytosis and environmental persistence by facilitating biofilm formation [McDougald et al., 2006; Williams et al., 2014a, 2014b]. We investigated expression of two genes involved in CPS production, wzb and cpsB, which encode a cognate phosphatase and mannose‐1‐phosphate guanyltransferase, respectively [Güvener and McCarter, 2003; Wright et al., 1990]. We also analyzed expression of a polysaccharide export periplasmic protein, SypC, which is homologous to a gene for EPS production in the “syp locus” of V. fischerii [Yip et al., 2005]. Twelve hours after entry into the estuary, wzb was downregulated, whereas both sypC and cpsB were highly upregulated (Figure 2a). In fact, cpsB was below the level of detection at t = 0 (indicated by the arrow on the graph) with the level of expression set at a 100‐fold increase. Generally, there was a decrease in expression of all three genes in response to changing oxygen conditions; however, this was not consistent through our sampling intervals.

Capsule production in the human host is an absolute virulence determinant, as cells lacking CPS succumb to the bactericidal effects of human serum, but its role in the environment is unknown [Williams et al., 2014b]. As production of these polysaccharides is metabolically expensive, there must be a key ecological role for their production in the estuarine environment, especially under hypoxic stress. Consistent with our previous in vitro results, the expression of wzb was decreased during conditions where oxygen was limited; however, the second CPS type and EPS production did not follow that trend [Phippen and Oliver, 2015b].

Having the ability to inhabit dynamic niches in a planktonic or a sessile state is an important factor for bacterial survival in estuarine environments. These lifestyles are driven by numerous factors, including but not limited to nutrient status, environmental conditions, and chemical signals. We found that the expression of pilA, a type IV pilin protein subunit which is required for type IV pilus formation [Paranjpye et al., 2007], was not driven by cyclical hypoxia in situ (Figure 2b). There was no significant expression change until 36 h in situ when its expression was decreased, followed by a slight increase, and the lowest level of expression was seen at 72 h (Figure 2b). Persistence of Vibrio species in the environment has been linked to the formation of another type IV pilin, mannitol‐sensitive hemagglutinin [Chattopadhyay et al., 2009], which we found to be upregulated up to 24 h and subsequently downregulated for the duration of our study (Figure 2b). Following the same pattern, flgE, coding for the monomeric hook subunit of the flagellum [Kim and Rhee, 2003], was also not influenced by hypoxic intervals (Figure 2b). Again, the decreased expression of both the type IV pili as well as flagellar synthesis after 24 h might suggest that V. vulnificus has already committed to an attached phenotype, which would also explain the initial increase in expression of mshA and flgE.

Quorum Sensing and Stress Response

V. vulnificus utilizes density‐dependent interspecific quorum sensing (QS) as a mechanism for impacting downstream regulation of genes involved in biofilm formation/dispersal, virulence, and in vivo survival [Miller and Bassler, 2001]. Additionally, QS has been shown to be integral in the dual life cycle of V. cholerae, driving both environmental survival and host invasion lifestyles [Zhu and Mekalanos, 2003]. We investigated the autoinducer‐2 (AI‐2) system of QS in V. vulnificus, which is currently the only known functioning system in this species [Liu et al., 2013], by measuring expression of the luxS gene responsible for production of the signaling molecule (AI‐2). We found that its expression was increased during hypoxic sampling events (Figure 2c) with over a tenfold increase at the first hypoxic time point. We also found that expression of rpoS, the alternate sigma factor involved in general stress response [Nowakowska and Oliver, 2013], was increased during hypoxia (Figure 2c). Furthermore, relA, responsible for the production of (p)ppGpp alarmones [Nowakowska and Oliver, 2013], was also upregulated in a cyclical hypoxic manner. It has been proposed that the expression of relA leads to the accumulation of RpoS, leading to increased protection from subsequent exposure to stress [Magnusson et al., 2005]. What is especially interesting was our observation that relA was below the limit of detection at t = 0 and then highly expressed at the initial hypoxic exposure, which could be why we did not see a robust change in rpoS at later time points (Figure 2c). The increased expression of multiple gene regulatory cascades suggests that hypoxia requires global and coordinated expression of many downstream products, which occur at regular intervals.

Virulence Factor Production

V. vulnificus has a large repertoire of putative virulence factors that contribute to survival in the host and subsequent colonization [Strom and Paranjpye, 2000]. Homologous to the multifunctional autoprocessing RTX (MARTX) toxin of V. cholerae, the RtxA1 toxin of V. vulnificus has a different mode of action, causing lysis of host cell membranes [Lee et al., 2007]. This toxin is transported out of the cell by a set of three genes (rtxB‐rtxD‐rtxE), and we utilized rtxB as a proxy for effective transport and function of the toxin [Lee et al., 2007]. Another regulator of virulence, the toxRS system, has been shown to increase the production of another putative virulence factor, vvhA, which encodes a potent hemolysin [Lee et al., 2000]. It has also been shown to regulate additional key virulence factors in V. cholerae and V. parahaemolyticus, and we examined toxR in this study [Lee et al., 2000; Lin et al., 1993]. Finally, we analyzed expression of a gene encoding a putative metalloprotease (metp), which our lab has found to be differentially expressed in human serum when compared to ASW [Williams et al., 2014c]. All three genes were found to be cyclically expressed (Figure 2d), with the highest differences occurring again at the first low‐DO sampling event. We also saw a large decrease (>100‐fold) in metalloprotease production at 48 h in situ, and overall, the expression of toxR varied significantly between low‐and high‐DO concentrations, whereas both rtxB and metp did not show significant changes at all intervals (Figure 2d). Although these factors are typically investigated in host pathogenesis, the high expression at times in situ, especially during hypoxia, leads us to hypothesize that they have distinct ecological functions.

Stressosome and Metabolism

We then examined the expression of the newly characterized stressosome module in V. vulnificus and metabolic gene representatives for both anaerobic and aerobic pathways. Recent analysis showed that V. vulnificus possesses genes homologous to the core of the stressosome, rsbR/S/T, which has generally been described as a gram‐positive mechanism [Pané‐Farré et al., 2005]. In Bacillus subtilis the stressosome senses environmental stress and transmits it to the alternate sigma factor σ B, which in turn regulates the expression of nearly 150 gene targets [Chen et al., 2003]. However, σ B is not present in the genome of V. vulnificus [Williams et al., 2014c]. Expression of all three genes, rsbR/S/T, was significantly upregulated in response to environmental hypoxia, as shown in Figure 2e. At each interval, there was differential expression upward of tenfold, suggesting that this module is influenced by environmental oxygen availability. Recently, Jia et al. [2016] described that the stressosome module in V. brasiliensis is activated by oxygen depletion and becomes inactive when oxygen is replenished, which is driven to be ligand‐binding dynamics [Jia et al., 2016]. The current study, in combination with the findings by Williams et al. [2014c] showing increased expression in ASW compared to human serum and by Jia et al. [2016], further validates the role of oxygen in the activation of the stressosome module in V. vulnificus. Lastly, we found that genes involved in both anaerobic (fnr) and aerobic (ilvC and purH) metabolism were upregulated in hypoxia (Figure 2f). The upregulation of fnr in low DO was anticipated as it is a main regulatory protein for the recognition of changing O2 levels and controls regulation of anaerobiosis [Spiro, 1994]. In addition, although both ilvC and purH are primarily used during aerobic respiration [Kim et al., 2003], they can also be utilized during anaerobiosis, which may explain their upregulated expression during hypoxia.

Correlation of Environmental Parameters and Gene Expression

To directly understand the relationship between environmental parameters, primarily DO, and expression of various genes in situ, we performed Pearson's correlational analysis on all expression results. A heat map of these correlations (Figure 3) revealed that there was a significant (p < 0.05) negative correlation between the expression of many of the genes investigated with respect to changing DO. Although we found no significant correlation of those genes involved in type IV pilin formation, flagella, CPS production (cpsB), alarmone synthesis, or two of the metabolic genes (ilvC and purH), they generally trended in a negative direction. Temperature (Figure 3, column 2) was significantly correlated with flgE, fnr, wzb, sypC, and two of the three stressosome genes; however, these genes were all positively correlated, as temperature was negatively correlated with DO. Both salinity and pH displayed very few significant correlations, except for toxR (pH and salinity), rsbR (salinity), and wzb (pH), all three of which were negatively correlated (Figure 3, columns 3 and 4). These results suggest that the primary driving condition for gene expression change was DO, although temperature did show many significant correlations.

Figure 3

Heat map of correlations between environmental conditions and gene expression by V. vulnificus. The correlation coefficients are indicated by color and range from −1 (red, negatively correlated) to 1 (blue, positively correlated) with 0 (purple) indicating no correlation between variables. Environmental values used were taken within ±1.5 h of the high‐ and low‐DO sampling times and were averaged at each sampling event. Averages of 2 biological replicates with 3 technical replicates each were used for expression values. Asterisks represent significant correlations (p < 0.05) between environmental parameters and gene expression (Pearson correlation).

Bacterial Community Composition

We collected water samples at each high‐ and low‐DO sampling event in order to assess bacterial communities, resulting in three high‐DO profiles and three low‐DO profiles. Figure 4a shows the relative abundance of the top 20 class‐level taxonomic categories identified, with the top seven categories representing approximately 75% of the total OTUs identified. Our analysis identified between 53 and 58 categories of classified OTUs at the class level, and those unclassified comprised less than 3.5% of all sequence data.

Figure 4

(a) Relative abundance of 16S rRNA sample libraries, (b) Bray‐Curtis measure of dissimilarity, and (c and d) beta diversity PCoA plots of bacterial communities collected at Hoop Pole Creek over a 3 day sampling event. Abundances represent the top 20 OTUs present in the samples. Bray‐Curtis represents relative abundances of all OTUs at one time point compared to the next time point; 0 = perfect similarity and 1 = complete dissimilarity. Beta diversity of both weighted (Figure 4c) and unweighted (Figure 4d) UniFrac distances is shown with colors corresponding to sampling date and DO profile.

In all samples, Synechococcophycideae accounted for approximately 26% ± 5% of all sequences, with the first sample (18 August 2015 low DO) having the lowest abundance at 16.7%, followed by Alphaproteobacteria (14% ± 1%), Oscillatoriophycideae (13% ± 3%), and Flavobacteriia (10% ± 1%). The high abundance of Synechococcophycideae is not surprising as it is a highly abundant picocyanobacteria inhabiting estuaries worldwide [Tan et al., 2015]. This group is responsible for a significant portion of estuarine primary productivity, and its abundance could in part explain why there was such a fluctuation in oxygen availability at this site [Murrell and Lores, 2004]. Additionally, previous studies have shown its abundance to be regulated most frequently by temperature and salinity, both of which remained relatively constant during our study [Sellner et al., 1988]. Furthermore, the top seven classifications are all regular estuarine inhabitants, which all fall under either subclasses of cyanobacteria, algal‐associated bacteria, or those commonly found in eutrophic sites.

We then analyzed β diversity from one sampling event to the next over the cyclical hypoxic intervals. For each time point, we calculated community dissimilarity between that sampling event and the next using the Bray‐Curtis dissimilarity measure, with 1 = perfect dissimilarity and 0 = perfect similarity. We found that the communities were highly similar, both temporally and in response to changing oxygen conditions, with only the first time point compared to the second being the quite dissimilar (Figure 4b).

Finally, we compared the community structures as a measure of both the relatedness of members in the community and what members were present. We calculated this based on weighted UniFrac β diversity measures after jackknifing the bacterial sequences. Visualization of the resultant principle coordinate analysis (PCoA) was performed in Emperor [Vázquez‐Baeza et al., 2013]. The weighted (Figure 4c) UniFrac analyses showed grouping of samples at low‐oxygen (purple and orange) and high‐oxygen (yellow and blue) concentrations, although there was one outlier of each group (red and green). It was interesting that the unweighted analysis (Figure 4d) showed low‐DO samples clustering together, whereas the high‐DO samples did not cluster with each other, possibly due to the sensitivity of rare sample abundance in unweighted analyses. This leads us to hypothesize that during normoxia there is a diverse community comprising many species, including those present at low abundance, and that diversity is subsequently lost as oxygen becomes limited. This lack of diversity during hypoxia should be further analyzed by deeper sequencing methods in order to better understand the specific differences between these conditions.

Conclusions

To our knowledge, this is the only comprehensive study on the role of in situ hypoxia on gene expression profiles of any Vibrio species. We observed the largely cyclical nature of these profiles to be negatively correlated with DO at this site (Figures 2 and 3). We also examined bacterial community compositions at both high‐ and low‐DO sampling events, finding that they were similar, although our unweighted β diversity suggested that rare sample occurrences might result in more diverse communities during normoxia (Figure 4). Importantly, this site experiences frequent nutrient inputs from a storm drain and represents a eutrophic estuarine environment, which should be explored further.

Coastal boundary systems and semienclosed seas include the Gulf of Mexico and the Baltic Sea, respectively, and are two areas which are likely to experience increased hypoxia due to global climate change projections [Barros et al., 2014]. The Intergovernmental Panel on Climate Change (IPCC) (2014) predicts a 4.37°C increase in the surface waters of the Baltic Sea from 2010 to 2099 and has already increased 1.35°C since 1980, whereas the Gulf of Mexico has only seen a 0.31°C increase from 1982 to 2006. This warming, combined with natural and anthropogenic eutrophication, increased microbial activity, and increased thermal stratification will likely exacerbate chronic and seasonal hypoxic events. Although hypoxia will increase globally via various mechanisms related to climate change [Barros et al., 2014; Rabalais et al., 2002, 2009], these two areas are interesting in terms of Vibrio infection rates, both historically and continuously from the Gulf of Mexico and newly emerging infections in the Baltic region [Baker‐Austin et al., 2013; Bier et al., 2015; Levy, 2015]. What is uncertain is the impact on diseases caused by marine Vibrios, in fact, the IPCC has low confidence in predicting the impact that this hypoxia has on these ecosystems, especially in the context of human health [Barros et al., 2014]. Although the National Oceanic and Atmospheric Administration (NOAA) and the International Council for the Exploration of the Sea monitor hypoxia in these areas, it has not been effectively applied to the increased incidence of Vibrio infections and should be included in future risk assessments.

These highly dynamic environments house a variety of opportunistic pathogens, including V. vulnificus, V. parahaemolyticus, and V. cholerae, all whose native ranges are expanding [Vezzulli et al., 2013]. This study highlights the importance of studying not only temperature and salinity as it relates to the abundance and ecophysiology of Vibrio species in situ, but that of oxygen availability as hypoxia will only increase as global temperatures continue to rise.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Supporting Information S1

Click here for additional data file.