Fluctuation and diversity of Hydromedusae (Hydrozoa, Cnidaria) in a highly productive region of the Gulf of Mexico inferred from high frequency plankton sampling

Medusa collection

Planktonic samples were collected, using a 100-micron net, 90 cm long, and a 30-cm mouth. The small mesh size prevented damage to any delicate hydromedusae that were collected. Plankton tows were conducted within the boat basin at Texas A&M University Galveston on Pelican Island (29°18′47.0′N 94°48′59.8′W). The basin receives seawater from Galveston Bay through the ship channel. Two tows per day were conducted three to four times a week from September 2015 to September 2016. The samples were collected during the morning, and each tow consisted of towing the net six times along the side of the dock by walking back and forth at a constant rate for a total of 156 m. The net was kept completely submerged in the water during the tows, to ensure that approximately the same amount of water was filtered for each sample. The plankton collected during the two consecutive tows were combined and considered as a single daily sample, and examined in the laboratory under a Leica M80 Stereomicroscope. Hydromedusae were isolated from other planktonic organisms using a pipette, counted and photographed using a Leica M80 Stereomicroscope connected to a Leica MC170 HD camera. Daily counts of medusae are shown in Fig. 1 and Table S1. Medusae were morphologically identified to the lowest possible taxonomic level using appropriate taxonomic keys (e.g., Bouillon et al., 2006) and preserved in ethanol for molecular analysis. The number of species and total hydromedusae abundance (i.e., individual count) of each sample was recorded.

Total medusa abundance and correlation with abiotic factors

Daily temperature (^∘C), salinity (ppt), dissolved oxygen (mg/L) (DO), and chlorophyll a (µg/L) of the Galveston Bay were made available by the Phytoplankton Dynamics Laboratory at Texas A&M University at Galveston. Temperature, salinity and DO were measured each morning at the same time as the plankton tows occurred; therefore, they reflect the water conditions at the time of sampling. Chlorophyll a data is available from January 2016 to the end of the sampling period only, so statistical analysis involving chlorophyll and medusa abundance was conducted on the January 2016 to September 2016 subsample (raw data are in Table S1). Descriptive statistics of the distribution of medusa abundance and the explanatory variables are reported for both the full sample (Table 1) and four subsamples based on the quartiles of the recorded temperature (Table 2). Medusa daily abundance (number of medusa collected in each sampling day) was plotted against the date of collection. A “bloom” was defined as any day with an abundance at least 1 standard deviation from the mean daily abundance (Miglietta, Rossi & Collin, 2008).

To investigate the potential non-linear effect of temperature on jellyfish abundance, we generated a bar chart of the total medusa count by temperature quartiles and a boxplot of the distribution of the abundance for the same temperature quartiles group (see Figs. 2A and Figs. 2B). The univariate analysis linking abundance to potential explanatory variables is concluded with the calculation of their correlation matrix (p-values testing the null hypothesis of zero-collation reported in parenthesis) (Table 3).

In addition to performing univariate analysis, we performed multivariate regression analysis based on ordinary least squares (OLS) in SAS. This approach allows for the identification of the effect of any of the explanatory variables on abundance, while controlling for the effect of the other explanatory variables included in the model. The general specification of the regression models is given by (1)${\displaystyle y_t={\alpha }_t+{\beta }^{{}^{\prime }}{x}_t+{\epsilon }_t}$ where t represent time (calendar day), y_t is abundance at date t, which is measured either in actual counts of medusa (model 1–5) or their natural logarithm (model 6–10), and x_t is a vector of explanatory variables (regressors) linked to medusa abundance (models differ depending on which regressors are included in the vector x_t). The use of the natural logarithm of medusa counts is proposed to account for the skewness of the distribution of daily counts and its potential detrimental effect on the statistical fit of the model. The elements of the vector β, which we estimate with OLS, capture the effect that the corresponding elements of the regressor vector (x_t) have on abundance. α_t is an intercept that is either constant (equal for every t) or varying by quarter (i.e., quarterly time fixed effects) depending on whether the model controls for potential unobserved factors not captured by the proposed explanatory variables. Finally, ε_t represent an error term uncorrelated with the regressors (x_t) that captures the unexplained variation in medusa abundance (Table 4).

Molecular barcoding

Species identification was confirmed using the hydrozoan barcoding molecule (a ∼600 bp fragment of the large ribosomal subunit of the mitochondrial RNA (lsu-rRNA, 16S)). Genomic DNA was extracted using a protocol modified from Zietara et al. (2000). The lsu-rRNA 16S was amplified using primers SHA (5′ACGGAATGAACTCAAATCATG T-3′) and SHB (5′-TCGACTGTTTACCAAAAACA TA-3′) and the following PCR conditions: 1 min at 94 °C, 35 cycles of 94 °C for 15 s, 50 °C for 1:30 min and 72 °C for 2:30 min, and a final extension at 72 °C for 5 min. PCR products were purified using exoSAP-it following manufacturer protocol. The purified PCR product was run on a 1% agarose gel stained with Sybrsafe at 100 W for 20 min to determine presence/absence of DNA. Confirmed PCR products were sent to the Genomics Core Lab at Texas A&M University Corpus Christi for sequencing analysis.

All sequences were edited in Geneious 10.0.5. Sequences from each species were run through the National Center for Biotechnology Information (NCBI) Basic Local Alignment Search Tool (BLAST) to confirm species identification. For each sequence and its most significant BLAST hit, identity scores and e-values were evaluated (see Table S2). Morphological analyses and barcoding and BLAST results were compared for correct species identification.

Species richness and dominant species

Species diversity was calculated monthly using the Shannon Weaver index. The species dominance index was calculated using Eq. (1) (Wang et al., 2016): (2)${\displaystyle Y=\frac{n_1}{N}{f}_i}$ where n is the number of individual species i; f is frequency of species i throughout the sampling period; N is the total number of individuals. Species with a dominance index more than 0.02 were taken as dominant species. The seasonality of five dominant species (Liriope tetraphylla, Blackfordia virginica, Malagazzia carolinae, Clytia gracilis, and Nemopsis bachei) and the dominant genus Obelia were further investigated using multivariate regression analyses. SAS University Edition was used to generate regression models to determine correlation with abiotic factors.

Time variation of total hydromedusa abundance

A total of 1,321 individual medusae were collected over 190 sampling days from September 1st, 2015 to September 30th, 2016 (Table 5). Samples were collected an average of 14.7 days per month. The daily average abundance for the sampling period was 6 medusae/day, the median 4, and standard deviation 11.3 (Table 1). The fact that the median is substantially lower than the mean suggests that the distribution of medusa counts is right-skewed, while the fact that the standard deviation is almost twice than the mean suggests that the variation is very high. Figure 3 reports the distribution of medusa counts and confirms the existence of infrequent days with extremely large counts. For instance, the maximum number of medusae found in a day was 104 and occurred on April 11, 2016. Over the 13 sampling months 19 blooms (defined as a day with an abundance at least 1 standard deviation from the mean daily abundance and represented in Fig. 1A by a dashed line) were recorded. The minimum abundance for a bloom was 19 medusae. Five blooms occurred during the summer, and were dominated by species such as Obelia, Clytia, Malagazzia carolinae, Liriope tetraphila. Six blooms occurred during the winter, and were dominated by Liriope tetraphilla, Nemopsis bachei and Ectopleura dumortierii. Eight blooms occurred during the spring and were dominated by Nemopsis bachei, Liriope tetraphilla (see Table 2 supplementary material for individual barcoded during each day). 49 out of 190 sampling days showed 0 medusa. The months with the highest medusa abundance were April (248 medusae), February (218), March (197) and September (182 medusae in 2015, and 173 in 2016). The months with the lowest medusa abundance were November (30 medusae), December (1), June (38), July (37) (see Table 5 for details). The presence of blooms, making the distribution of medusa counts right-skewed, has ultimately prompted us to propose several regression specifications with the dependent variable expressed in logarithms.

Hydromedusa and environmental factors

Temperature varied from 11.7 °C (min) to 32  °C (max), salinity from 6.00 ppt to 31.97 ppt, and DO from 4.23 µg/L to 10.42. µg/L , chlorophyll a from 0 to 9 µg/L (Table 1). Data for chlorophyll a (µg/L) are available each sampling day from January 2016-September 2016 only (Fig. 1C). As can be seen from Table 3, no significant correlation was found between the number of medusae and temperature. However, inspection of the data suggests that the relation between temperature and abundance is not linear, and a simple correlation measure would fail to capture such non-linearity. Figure 2 reveals that medusa are abundant during periods of mild temperature. Specifically, the bar chart shows that, by dividing the sample into 4 equally-sized temperature quartiles, more medusae were captured during days of intermediate temperatures (the two middle quartiles) relative to days with more extreme weather (first and fourth quartiles), the difference being in the order of a 50% increase. The boxplot in Fig. 2B further explores the distribution of medusa abundance over the same temperature sub-periods and confirms that the median counts (horizontal bar within the boxes, which capture the mid-50% of the distribution) are also higher during the days of mild temperature. Moreover, the boxplots reveal that positive outliers (typically blooms) are more likely during mild temperature days.

The correlation between the other explanatory variables and abundance also seems weak (Table 3), with salinity being marginally statistically significant and remaining variables being insignificant. However, a univariate analysis linking environmental factors to abundance individually is likely deficient because these factors likely interact to produce their outcome on medusa abundance. Moreover, given the high correlation among the explanatory variables, univariate analysis could result in an omitted-variable bias.

To better understand the effect of environmental factors on medusa abundance, multivariate ordinary least squares (OLS) regressions were run using the daily medusa abundance (both in levels and natural logarithms) as the dependent variable, and temperature, salinity, DO, and chlorophyll a as the independent variables. The first 5 models use actual medusa counts as the dependent variable, while the next 5 models use the natural logarithm of counts. To account for potential non-linear effects of temperature on abundance, we either include a quadratic term (Models 2–4, 7–9) or several temperature quartile dummies. As can be seen from the table, models ignoring non-linearity (Model 1 and 6) have a relatively poor model fit, as evidenced by a lower R squared. Moreover, these naive models would lead to the erroneous conclusion that temperature is unrelated to abundance due to lack of statistical significance (high p-values). On the other hand, the other models in the table point to a strong and statistically significant non-linear effect of temperature on abundance. For instance, the quadratic model (column 2 in the table) shows that both the linear and quadric term are statistically significant at the 5% level. The interpretation of the coefficient is that for low levels of temperature the linear term (i.e., 2.16) dominates, but the quadratic term (−0.08) becomes more important as temperature rises. These estimates indicate an intermediate level of temperature associated with maximum abundance. Specifically, this level is obtained by solving 2.16–2*0.08*Temperature = 0, which is approximately 21 °C in this case. A somewhat more general approach, one that does not rely on a stringent quadratic assumption, for assessing the impact of temperature is to use several temperature dummies. The estimates in column 3 suggest that, relative to the benchmark case of high temperature, mildly lower and mildly higher temperatures are associated with respectively 8.7 and 6.5 additional individuals per sampling day on average. Considering the fact that on average 5 individuals are collected during hot days (not reported), our estimates represent an average increase of more than 100%.

To account for potential seasonality and unobserved factors driving medusa abundance we also run some regressions (Model 3 and 8) including quarterly time fixed effects (time-varying intercept, α_t, in Eq. (1)). The estimates of time fixed effects, which we do no report in the table for brevity, are obtained by including calendar quarter dummies in the regression. As can be seen (Model 3 and 7), the inclusion of time fixed effects increases model fit (higher R squared), but does not drive away the impact of the other regressors (especially temperature) on abundance. This result is important because it helps us rule out with some confidence that the impact of temperature we document is not spurious, or driven entirely by unobserved factors that we might not be accounting for. Finally, while the coefficient estimates on temperature and squared temperature in Model 3 appear to differ from those of Model 2, we note that a similar reasoning used above reveals an approximate level of temperature of 20.5 degrees Celsius associated with maximum abundance. Thus, our finding on temperature are robust across different models.

Species identification and molecular barcoding

Morphological identification of planktonic hydromedusae is difficult because of their small size, the high number of species, and the morphological disparity through out the life cycle, as new-born medusae and adult medusae may look substantially different. Identifying hundreds of medusae to the species level becomes a time and cost-prohibitive enterprise. To manage species identification we thus adopted a combined approach between morphological identification and molecular barcoding.

After collection and during sorting, medusae (a total of 1,321) were identified at the lower possible level using appropriate keys. Most medusae could be identified to the genus level, some at the family level, only few could be identified at the species level. Specimens were also photographed, and fixed in ethanol for barcoding studies. DNA was then extracted from all the individual medusae and 470 mitochondrial ∼600 bp sequences of the large ribosomal subunit of the mitochondrial RNA (lsu-rRNA, 16S) were generated (See Table S2, for Blast results and Genbank accession number of the newly generated sequences). This represents 36% of all specimens. The relatively low success rate is possibly due to the small size of some of the medusae (i.e., Obelia). Each sequence was compared with the online repository Genbank using the BLAST feature and results were used to augment morphological identification. For each query sequence, blastn of the NCBI-BLAST package was run locally, and the result was parsed to retrieve topmost match. We grouped the sequences in 3 categories depending on the level of identity with their best match as follow: identification through Genbank was considered reliable when genetic identity with best match was between 100% to 98% (category 1); it was considered acceptable but with lower confidence, when genetic identity was between 97.9 to 95% (category 2); identification was considered uncertain when genetic identity was lower than 94.9% (category 3). These specimens probably belong to the genus or family of best match. We understand these are arbitrary threshold and need to be taken as an aiding tool to morphological identification. They also can help spotlight the amount of species that are currently missing in Genbank.

Table S2 shows the 470 individual medusa that were successfully sequenced and their BLAST results (best match, % identity between sequences, e-value and query coverage). Using the combined approach between morphological analyses and barcoding we identified 967 individual medusae as belonging to 25 species (see Table 5 for a complete list). Of this 11 species, or 44% of total (Aequorea victoria, Blackfordia virginica, Clytia folleata, Clytia gracilis, Malagazzia carolinae, Nemopsis bachei, Obelia bidentata, Obelia dichotoma, Turritopsis dohrnii and Sertularella diaphana) were in category 1 (identified with confidence), 6 (24%) (Bougainvillia triestina, Bougainvillia muscus, Stauridiosarsia reesi, Ectopleura dumortierii, Hydractinia americana, Liriope tetraphylla) in category 2 (identified with lower confidence), and 8 (32%) in category 3 (uncertain identification). These were mostly rare species found only once and with relatively low number of individuals (their best match in Genbank was Clytia elsaeoswaldae, Corymorpha nutans, Clytia sp. 1, Earleria quadrata, Eucheilota maculata, Koellikerina fasciculata, Lovenella assimilis, Obelia geniculata). 354 individual medusae for which barcoding sequence was not available and that could not be identified morphologically (because undescriptive juvenile or damaged medusa), remain unidentified (see Table 5). Sequences are in Genbank under the following accession numbers: MN364418 –MN364461, MN355950 –MN355952, MN355770 –MN355845, MN355846 –MN355904, MN355575 –MN355691, MN355905 –MN355911, MN355912 –MN355949, MN355718 –MN355768, MN355692 –MN355717, MN355508 –MN355516, MN347011 –MN347014, MN343754 –MN343758, MN341844 –MN341854 and MN367046 –MN367055.

The community composition and seasonality of dominant species of hydromedusa

Malagazzia carolinae was the most common species, found in 10 out of the 13 sampled months. Several species including Aequorea australis, Bougainvillia muscus, Clytia elsaeoswaldae, and Turritopsis dohrnii were found only once during the sampling period. Many Obelia medusae for which barcoding was not successful, were morphologically identified to the genus level (due to their distinctive generic characters) but could not be identified at the species level, and were all collectively called Obelia spp. It is possible that this category contains more than one Obelia species. Table 6 shows the dominant species for each month, and Table 7 shows the diversity index (Shannon Weaver) per month. Only one medusa was found during December 2015 and could not be identified to species, so it is not included in the dominance analysis. The genus Obelia dominated the months September 2015–November 2015 and August 2016–September 2016. Nemopsis bachei was the dominant species January 2016–March 2016, Liriope tetraphylla was the dominant species in April 2016. The dominant species in May 2016 was Blackfordia virginica. Malagazzia carolinae dominated June 2016–July 2016. March, July and September were the months with higher diversity index (Table 7). March 2016 had the highest number of species (14) followed by July and September 2016 (11 each). December 2015 and January 2016 had only one species (Fig. 4 and Table 5).

Relationship between dominant species and environmental factors

Figure 5 show the plots of each species abundance with temperature, salinity, and DO and suggests that each species has a distinct seasonality. Regression models were run for the top five dominant species (Blackfordia virginica, Clytia gracilis, Liriope tetraphylla, Malagazzia carolinae, Nemopsis bachei), and the genus Obelia, to determine their specific relationships with environmental factors. The species-specific models were run using medusa abundance as the dependent variable, temperature (linear and non-linear), salinity, and dissolved oxygen as the independent variables, and time fixed effects classified by season (Table 8). This corresponds to Model 3 in total medusa count (see above and Table 4). Model 3 was chosen because it incorporates linear and non-linear temperature regression and time fixed effect. Under this model, temperature had a significant effect on the abundance of Clytia gracilis, Liriope tetraphylla, Malagazzia carolinae, and Obelia. Salinity and DO, but not temperature, had a significant effect on the abundance of Nemopsis bachei.

Environmental variables

Temperature has a non-linear relationship with medusa abundance (Table 4). Using regression analysis based on a quadratic model (Table 4), we have identified a level of temperature of approximately 21 °C that is associated with maximum medusa abundance in Galveston Bay. Regression analysis using temperature quartiles as dummies produces finding consistent with the quadratic model: the highest total medusa abundance is found during days with moderate temperature (medium low and medium high temperature), and decreases sharply otherwise (Figs. 2A and 2B). Temperature also shows a strong effect on the abundance of many of the dominant species (Table 8). This is consistent with previous studies conducted on individual species of Hydrozoa (Ma & Purcell, 2005; Nowaczyk et al., 2016; Wintzer, Meek & Moyle, 2013). Salinity was not significantly correlated with total medusa abundance, however Nemopsis bachei showed seasonality and had a significant relationship with salinity and DO (Fig. 5). This is consistent with previous studies on this species, which have shown a correlation with salinity and not temperature (Nowaczyk et al., 2016). Blackfordia virginica was also one of the dominant species during spring but was not statistically correlated with any of the tested environmental factors (Table 8). These individual species analyses and the fact that dominant monthly species changed with time (Table 6) indicate that trigger of medusa production may be species-specific (Nowaczyk et al., 2016; Wintzer, Meek & Moyle, 2013).

Gulf of Mexico and Galveston Bay Hydrozoa diversity

The most recent species list of Hydrozoa for the Gulf of Mexico was compiled by Calder & Cairns (2009) and listed 214 species. Only eight of the species recorded in this study were on the Calder & Cairns checklist. An additional seven species had been reported in the GoM in previous studies (See Table 9). In total, eight species found in this study had never been recorded in the GoM before. The most recent study on Hydrozoa in Galveston Bay was conducted by Defenbaugh and Hopkins in 1973. They surveyed only polyps and found 25 species. Only 4 of the 25 species found in this study were previously described in Defenbaugh and Hopkins (Table 9). Thus 21 species are recorded for the first time in the Galveston Bay.

Ours study focuses on the medusa stage only, thus neglecting Hydrozoa that have lost the medusa stage (with the exception of one species Sertularella diaphana, that was caught in our plankton tows in the form of sporosacs). This may explain some of the discrepancies between the species found in our sampling and those found in other published studies that focused only on benthic sampling. It is also noteworthy that the Galveston Bay area has a large amount of ship traffic exposing it to potential species introduction through ballast water, which may dramatically alter its species diversity and composition through time (Steichen et al., 2012). The majority of the dominant species (with the exception of Malagazzia carolinae) found during this research are however widely distributed throughout the GoM (namely Blackfordia virginica, Clytia gracilis, Liriope tetraphylla, Nemopsis bakei) (see Table 9 for references). In contrast, the dominant Malagazzia carolinae was first described in 1900 by Mayer in Tortugas, Florida, an island at the edge of the Gulf of Mexico, and never recorded in the GoM. M. carolinae is generally found on the coasts of New Zealand and China (Bouillon, 1995; Du et al., 2011). More surveys will be required before we can know the extent of the distribution and persistence of Malagazzia carolinae in the GoM. Aequorea australis, Bougainvillia triestina, Stauridiosarsia reesi, are new to the GoM. The introduced Turritopsis dohrnii (Miglietta & Lessios, 2009; Miglietta, Maggioni & Matsumoto, 2018) was recently described for the first time in the GoM. All of these species had BLAST identity above 95% (Table S2), but had sample sizes lower than 5 individuals. Other species with low BLAST identity (thus with uncertain identification) were found for the first time in the GoM , and their BLAST top hits were Clytia elsaeoswaldae, Earleria quadrata, Eucheilota maculata, Koellikerina fasciculata, and Lovenella assimilis and Obelia spp. Those are rare and/or understudied species, as shown by the lack of high identity matches on Genbank. More data from the GoM and worldwide and in-depth morphological identification are needed in order to further discuss the implications of these findings. Our results however indicate that the Hydrozoa in this region of the world are still heavily underestimated. They also show that taxonomic studies that couple molecular barcoding and morphological identification, and studies that focus on the full life cycle of hydrozoan species (polyp and medusa) are needed to assess biodiversity. Monitoring local biodiversity and its long-term changes due to human disturbance is particularly important in high traffic marine areas such as Galveston Bay and the GoM in general.