Impacts of hypoxia on estuarine macroinvertebrate assemblages across a regional nutrient gradient

Because of its relatively low salinity, the Gulf of St. Lawrence has been referred to as a very large estuary (Therriault 1991), but this study focuses on the many small estuaries within it, specifically those situated in the southern portion (henceforth referred to as the southern Gulf of St. Lawrence, SGSL), with watershed areas ranging from 37 to 386 km² in this study (Fig. 1). Low freshwater input results in estuaries that are well mixed throughout the water column (Bugden et al. 2014) and that remain saline throughout. Tides vary in both amplitude and nature and can be semi-diurnal or diurnal (Pingree and Griffithis 1980; Godin 1987; Koutitonsky et al. 2004). Because of local geography and low freshwater input, estuaries are generally shallow and small. On the north shore of Prince Edward Island most estuaries are lagoon-type, with barrier islands, and drain into the southern Gulf, whereas those on the south shore empty into the Northumberland Strait and are most often coastal embayments (Glibert et al. 2010). Eelgrass (Z. marina), the local seagrass, is in decline in the region (DFO 2009), particularly in the upper reaches of Prince Edward Island (PEI) estuaries (Hitchcock et al. 2017). Many PEI estuaries have agriculturally driven high nitrate–N loads (Danielescu et al. 2007; Danielescu and MacQuarrie 2011), are dominated by sea lettuce (Ulva spp.), and experience seasonal anoxia as a result (Bugden et al. 2014). The remaining provinces in the southern Gulf region have ≤5% of their land use in annual crops (Natural Resources Canada 2009; Grizard 2013). Thus, the estuaries of the region constitute a gradient of nutrient enrichment through which to examine the interactions of factors structuring invertebrate communities (Table 1). The transition zone from 0 to 15 PSU salinity is short and creates a highly susceptible area to nutrient impact immediately downstream, and is also the area of interest in this study. This area is algae-dominated at nutrient-impacted sites (Bugden et al. 2014) and eelgrass-dominated at oligotrophic sites. At sites that were neither eutrophic nor oligotrophic (i.e., along the nutrient continuum), vegetation tended to be dominated by one or other type but was patchier. In total, 13 estuaries spanning a gradient of nutrient loading throughout the SGSL and that were dominated by either Ulva or Z. marina, were sampled monthly for epifauna and infauna and more frequently for standard physicochemical variables, e.g., salinity, temperature, and pH. Dissolved oxygen was the focus of this study and was continually measured using optical logger technology.

Water chemistry variables were measured biweekly, before noon and within 2 h of slack tide, 0.5 m from the substrate using a YSI V2 6600 multiparameter sonde (Yellow Springs Instruments, Yellow Springs, Ohio, USA). The sonde was equipped with optical probes for dissolved oxygen, pH, conductivity, and temperature. Environmental variables were taken biweekly from June to September 2013 and the overall averages are presented in Table 1. For the monthly invertebrate sampling period, however, these variables were averaged over the two biweekly samplings preceding it to represent typical conditions prior to sampling. One litre grab samples of water were taken for chlorophyll sampling concurrently with environmental variable measurements, about 0.5 m below the water surface, and stored in the dark until processing later that day. Water samples were filtered through glass filter fiber paper type F and then submerged in 5 mL of acetone to extract chlorophyll. Samples were stored in the dark at −80 °C and later analyzed using high-performance liquid chromatography (Schein et al. 2012).

An Onset HOBO (Bourne, Massachusetts, USA) dissolved oxygen logger, which recorded dissolved oxygen (mg/L) and temperature (°C) hourly, was used to capture fine-scale variability in the oxygen regime. Criteria for the area of deployment for the dissolved oxygen logger, and later for invertebrate sampling, were as follows: salinity (average PSU of 15–25), depth (∼1.5 m at slack tide), and vegetation presence (Z. marina or Ulva). Cumulative dissolved oxygen metrics for the 48 h preceding sampling were selected to capture symptoms of nutrient impact, i.e., hypoxia and dissolved oxygen supersaturation. Metrics were based on work by Coffin et al. (2018), which found that hypoxia (proportion of time below 4 mg/L), dissolved oxygen supersaturation (proportion of time above 10 mg/L), and the coefficient of variation of dissolved oxygen were all predicted by water residence time and nitrate loading. Additionally, metrics for the proportion of time <2 mg/L, <6 mg/L, and >15 mg/L were incorporated for greater resolution of sites that were minimally and maximally impacted by nutrients (Miller et al. 2002; Landman et al. 2005; Riedel et al. 2014; Coffin et al. 2018). The duration of 48 h was selected because effects from hypoxia prior to invertebrate sampling would likely still be evident (Coffin et al. 2017), and although not investigated directly in this manuscript, the impacts on behaviour and survivorship can occur after extended exposure to hypoxia (Miller et al. 2002; Riedel et al. 2014). There were two data quality issues over the course of the sampling in the Murray and Cocagne estuaries. Data from the Murray logger were unreliable as it was periodically buried in anoxic substrate and were, therefore, excluded from analyses or ordinations involving environmental variables. The logger from Cocagne, a Z. marina site, became fouled with vegetation immediately prior to the July sampling and so the first 48-h period without fouling was used instead (approximately 5 d earlier).

A tidal prism model employing estuarine volume, mean tidal amplitude, and freshwater input was used to create a proxy for water residence time. Bathymetric data were collected in a related project and used here to calculate estuarine volume using ArcGIS 10.3. Tidal data were collected every 10 min for at least 30 d at the dissolved oxygen logger location in each estuary using Onset level loggers (one barometric and one submerged). Harmonic tidal models were created from those data using the t_tides program in Matlab (Pawlowicz et al. 2002). These models were used to simulate tides over the period of dissolved oxygen logger deployment (15 May to 30 November 2013), and mean tidal amplitude was derived by averaging the tidal amplitudes that occurred over that period. Higher residence time values are indicative of sites with low flushing rates. These values represent the ratio of overall estuarine volume relative to freshwater input and tidal exchange.

Nitrogen loading (nitrate-N) data were collected from two sources: primarily from Grizard (2013) and supplemented with data from Bugden et al. (2014). In short, daily water flow (m³/s) was multiplied by mean nitrate-N concentration (kg/ha/year), which is the most bioavailable form of nitrogen and best correlated with plant productivity (Hemminga and Duarte 2000). A strong relationship was found between daily flow and watershed area for rivers in this region (n = 8, R² = 0.97). As our area of interest was farther downriver than the most downstream measurements, the calculated nutrient loads, and flow, were prorated for the larger watershed area and assumed to be proportional to the measured watershed area.

A Birt–Flannagan-modified Ekman dredge (15 cm × 15 cm) was used to collect sediment samples at three locations, but only on one occasion, within the sampling region for particle size and organic content analysis. Sediment was homogenized and then dried at 60 °C for 48 h. A 5 g sediment subsample from each replicate was analyzed for organic content by loss on ignition at 550 °C for 24 h. The organic-free sediment sample was analyzed for grain size using laser diffraction (Horiba laser particle sizer model LA-960). During the analysis, ultrasonic dispersion was used to ensure homogeneity of the sample. For our analyses, the 10th, 50th, and 90th percentiles were highly correlated within estuaries, and only the 50th percentile (D50) was retained to characterize grain size.

Invertebrate sampling occurred before noon and very close to slack tide, monthly between 10 June and 25 September 2013. Sampling was boat-based, in 1–2 m of water, and was within a 200 m radius of the dissolved oxygen logger (location as described above). To ensure the independence of the samples, ephemeral algal mats were sampled by traveling in a random direction, mooring, and then attempting collection. Several sampling methodologies were employed to suit the assemblage or habitat being sampled (Z. marina, Ulva, or infauna). Epibenthic invertebrates, henceforth referred to as epifauna, in Z. marina habitat were sampled using a Birt–Flannagan-modified Ekman dredge (15 cm × 15 cm) dropped into a Z. marina bed. This method was ineffective for Ulva habitat as the dredge was incapable of cutting through the thick algal mats. Thus, an alternative method was devised using two modified bow-head garden rakes (40.6 cm long with 2.5 cm between tines). The sampler would place the rakes, one in each hand and shoulder width apart, on the substrate and bring them together underwater, sampling an approximate area of 0.25 m² (Coffin et al. 2017). Vegetation samples for either habitat were brought into the boat and into a bucket of invertebrate-free water. Vegetation was visually inspected for degradation (only healthy algae were retained), separated from sediment and manually cleared of invertebrates, and then placed into a plastic bag to be processed in the lab. Ulva samples were large and compact when brought into the boat, and a total of 135 Ulva samples were taken with an average dry mass of 20.9 ± 1.5 g (SE); although some invertebrates on exposed edges may have escaped during sampling, the majority were contained within the algae ball. The remaining water–sediment–macroinvertebrate slurry was immediately sieved through 500 μm mesh and stored in 95% ethanol to be processed later. Infauna were collected using a 7.6 cm diameter core, about 15 cm deep within 5 m of the vegetation samples, also sieved through 500 μm mesh and stored in 95% ethanol. For each sampling method, five replicate samples were taken each month (June, July, and August/September). All samples were further processed using a dissecting microscope (40× magnification), and macroinvertebrates were identified using appropriate taxonomic guides (Bousfield 1973; Appy et al. 1980; Pollock 1998; Merritt et al. 2008; Thorp and Covich 2010). As sampling methodologies differed for invertebrate collection, all analyses were performed using relative abundances. A selection bias for different fauna was likely given the variety of sampling equipment used, and thus infauna and epifauna were not compared statistically. Vegetation was not measured explicitly in this study, but Z. marina coverage estimates were taken in a related study (Hitchcock et al. 2017). Within the sampling area, only Montague and Wilmot had sparse Z. marina near the sampling area at Ulva-dominated sites and none of the Z. marina sites had significant Ulva incursion (Hitchcock et al. 2017). Generally, the area sampled had continuous vegetation; sites where bare substrates existed were not sampled.

All data were analyzed using Plymouth Routines in Multivariate Ecological Research package version 6.1.18 (Clarke and Gorley 2006) with PERMANOVA+ (Anderson et al. 2008) and STATISTICA version 12. Environmental variables from the two weeks prior to each sampling time at each site were averaged and then analyzed using principal component analysis (PCA) to examine relative differences between sites. To meet assumptions of linearity and homoscedasticity for parametric analyses, environmental data were tested visually and with Cochran’s C test, respectively, and, where necessary, log transformed and then normalized for subsequent analyses (Clarke and Gorley 2006). For invertebrate assemblage data, Bray–Curtis similarity resemblance matrices were created with the inclusion of a dummy variable to account for samples with no species (Clarke and Gorley 2006; Anderson et al. 2008). Principal coordinate analysis (PCoA) ordinations were created, using the Bray–Curtis dissimilarity matrices, for epifaunal and infaunal communities to visualize data by habitat type. Mixed-model PERMANVOAs were used to test between the fixed factor habitat type (Z. marina or Ulva) and “Month” (the convention for repeated measures designs (Anderson et al. 2008)), and the random factor “Site” nested within “Habitat Type”. Additionally, species richness was calculated for each site and sampling time and then analyzed for each sampling methodology by habitat type using two one-way ANOVAs. Next, habitat types were analyzed independently using mixed-model PERMANOVAs to examine the effects of the fixed factor “Month”, the random factor “Site”, and the potential interaction between those factors. Sites were designated as a random factor because they were chosen from a larger group of potential estuaries on which inferences were to be made. In Z. marina habitat for both sampling methodologies, initial analyses resulted in a negative estimate for the component of variation and large p values (p = 0.542 and p = 0.479 for epifauna and infauna, respectively) for the temporal factor “Month”. As suggested by Underwood (1997) and Anderson et al. (2008), this term was pooled with the interaction term, “Month × Site” in this case. Similarity percentages were also calculated across all sites and sampling times for Z. marina and Ulva habitats to determine which species contributed most to the dissimilarity.

Relationships between invertebrate assemblages and environmental predictor variables were examined using distance-based linear models (DISTLM) and visualized using distance-based redundancy analysis (dbRDA). Correlations between independent variables were restricted to r < 0.7 (Pearson’s correlation), as strong correlations have undue influence over results and can lead to spurious interpretation of the dependent variables (Anderson et al. 2008). When independent variables were highly correlated, the variable that was most correlated with other independent variables was eliminated. DISTLM is used to analyze the relationship between the Bray–Curtis resemblance matrix and potential explanatory variables, in this case invertebrate assemblages and associated environmental variables. Because some of the environmental data were at the monthly scale, invertebrate data were pooled (averaged) by sampling time for each estuary to enable analyses. As with the mixed-model PERMANOVAs, analyses were conducted on each habitat type and for epifauna and infauna independently. Model selection for the DISTLM was performed using step-wise selection based on adjusted R² criteria.

Environmental variables from the 13 estuaries studied, except for Murray because of logger burial in anoxic substrate (Table 1), were visualized using a PCA (Fig. 2). The first three axes of the PCA explained 33.9%, 21.1%, and 19.8% of the total variation, respectively (Table 2). There was separation of sites that drain into the southern Gulf and those that drain into the Northumberland Strait. This can be explained by SGSL-draining sites tending to be more correlated with longer water residence time and higher values for variables associated with eutrophication, specifically chlorophyll, dissolved oxygen variability (DO CV), and hypoxic hours (dissolved oxygen <6) (Fig. 2). The primary axis was driven by variables related to nutrient impact (e.g., coefficient of DO CV, the percentage of hours of dissolved oxygen under 6 mg/L, organic content of the sediment sample, water residence time, and chlorophyll) (Table 2). Salinity was also generally higher at Ulva sites, most likely due to their dominance in SGSL-draining estuaries that have lower tidal variability compared with Northumberland Strait draining sites (Table 1). Conversely, the secondary axis was driven by nitrate-N loading, average substrate grain size (D50), and organic content (Fig. 2, Table 2). The tertiary axis was nearly as important as the secondary axis and, like the primary axis, was driven by eutrophic variables, specifically dissolved oxygen metrics related to high and low oxygen and also pH; temperature was not important for any of the first three axes (Table 2).

The PCoA ordinations examining relative abundance of epifauna and infauna between Z. marina and Ulva habitat showed similar patterns and in both cases the first two axes explained more than 50% of the total variation (Fig. 3). Samples were generally grouped according to vegetation type within the ordinations for epifauna and infauna. It is notable, however, that Kouchibouguac, a Z. marina-dominated site that had a high proportion of hydrobid snails, was more associated with Ulva sites. For both epifauna and infauna sampling methodologies, hydrobid snails dominated Ulva habitat but were also found regularly at low densities in Z. marina. Conversely, cerithid snails were uncommon in Ulva habitat but dominated in Z. marina (Table 3). These two taxa contributed most to the average dissimilarity between the invertebrate communities found in the two plant habitats for both epifauna and infauna. Epifauna were further distinguished by Gammarus mucronatus Say, 1818, Gammarus lawrencianus Bousfield, 1956, capittelid polychaetes, and the softshell clam Mya arenaria Linnaeus, 1758. Infauna were distinguished by M. arenaria, capittelid polychaetes, Acteocina canaliculata (Say, 1826), and nereid polychaetes (Table 3). Nested PERMANOVAs revealed that invertebrate communities were significantly different between vegetation types and sites but not sampling times (Table S1). There was a significant interaction between Site (vegetation type) × Month; a posteriori tests revealed that each site was significantly different from all other sites for every month (results not shown).

For both sampling methods, pooled across time and site, there was significantly higher species richness in Z. marina than Ulva habitat (ANOVAs for epifauna and infauna: F _1,37 = 20.182, p < 0.001 and F _1,35 = 7.744, p < 0.009, respectively; Table 4). The different sampling methodologies precluded a statistical comparison of abundances which were higher in Ulva habitat. Gammarid amphipod abundances, for example, differed by up to two orders of magnitude (1000s vs. 10s) in Ulva relative to Z. marina sites. Overall, hydrobid snails, several species of Littorina Ferussac, 1822, Tritia obsoleta (Say, 1822) (previously Ilyanassa obsoleta), G. mucronatus, G. lawrencianus, and corophid amphipods were more abundant in Ulva habitat, whereas cerithid snails, A. canaliculata, Astyris lunata (Say, 1826), and capitellid polychaetes were more abundant in Z. marina (Table 3).