Shoreline oiling effects and recovery of salt marsh macroinvertebrates from the Deepwater Horizon Oil Spill

Total Petroleum Hydrocarbon (TPH) analysis

Surface soils (0–2 cm) were collected from each station on each collection date. Samples were extracted with dichloromethane (DCM), and analyzed gravimetrically (Lin & Mendelssohn, 2012). DCM extracts were transferred to pre-weighed dishes, where the DCM was evaporated, and unevaporated oil remaining in the dishes was weighed to the nearest 0.0001 g. TPH concentration was calculated and expressed as mg g⁻¹ dry soil.

Aboveground and belowground plant biomass

Samples for aboveground and belowground plant biomass were taken within a haphazardly located 0.25 m² quadrat at each station. All plants rooted within the quadrat were clipped to the ground surface and separated into live and dead components by species (mostly S. alterniflora and J. roemerianus). All aboveground biomass was then dried to a constant mass at 60 °C and weighed. Belowground biomass was determined by extracting a 7.62-cm diameter by 36-cm long soil core from each station on each sampling event. Cores were subsequently sectioned into 6-cm segments and washed over a 2-mm mesh sieve to remove sediment and particulate organic matter. Material from each segment was then separated as live and dead roots and rhizomes, and live material was dried to a constant mass at 60 °C, and weighed. Here we present data from the uppermost 6 cm. See Lin & Mendelssohn (2012) and Lin et al. (2016) for additional details.

Macroinvertebrate sampling

To quantify fiddler crab (Uca) species composition, burrow density, and burrow diameter, and marsh periwinkle (Littoraria irrorata) density and shell length, three 0.25 m² quadrats were haphazardly placed approximately 1 m landward from the marsh edge at each station and within 5 m of the station location stake. We collected L. irrorata from the marsh surface or attached to the vegetation within the quadrat; juvenile L. irrorata were often discovered in the leaf bracts and furled, dead leaves of S. alterniflora. L. irrorata shell length, from the base of the aperture to the top point of the shell, was measured to the nearest millimeter using calipers. Once measured, L. irrorata were returned to the sampled quadrat.

L. irrorata were grouped into putative age classes based on shell length: snails <6 mm were classified as juveniles, snails 6–13 mm as sub-adults, and snails >13 mm as adults (Bingham, 1972; Hamilton, 1978; Stagg & Mendelssohn, 2012; Zengel et al., 2014; Zengel et al., 2015). We additionally classified snails >20 mm as “large adults”. Length-frequency distributions of L. irrorata are often bi- or tri-modal with the modes representing age class (Hamilton, 1978; Pennings et al., 2016; Zengel et al., 2014; Zengel et al., 2015). Based on well-established growth estimates (Stiven & Hunter, 1976), juveniles are considered to be <1 year of age, sub-adults ∼1 year, and adults 2 years and older (Pennings et al., 2016). Evaluating changes in shell size over time by comparing size frequency histograms has previously been used as an aid to determine oiling impacts and to estimate recovery trajectories (Pennings et al., 2016; Zengel et al., 2014; Zengel et al., 2015; Zengel et al., 2016c).

Because fiddler crabs are mobile and direct observation is difficult, we used burrow density as a proxy for population density. A strong relationship exists between fiddler crab carapace width and burrow diameter (Mouton & Felder, 1996), and burrow diameters were measured to the nearest millimeter with a transparent ruler (after Pennings et al., 2016). When present in the sampling site vicinity, fiddler crabs were captured, identified to species, and released. We were not able to count or measure crab burrows on four sampling occasions (36, 40, 60, and 66 mo after the spill) because the marsh platform was tidally inundated.

Statistical analysis

All statistical analyses were conducted using SAS (Statistical Analysis Systems, version 9.2, SAS Institute, Cary, NC). Dependent variables (Uca burrow density; Uca burrow diameter; total, adult, subadult, and juvenile L. irrorata density; and L. irrorata average shell length) were each tested with repeated measures Analysis of Variance (ANOVA) Mixed Models (with oiling category and time as the main effects and with an oiling category * time interaction), with pooled variance and with natural log transformation. Natural log transformed data with pooled variances met the assumption for normality (based on the Shapiro–Wilk test) for most tests. Tukey’s test was used to make pairwise comparisons.

Because qualitative trends associated with periwinkle length and density were noted, a quadratic regression mixed model Analysis of Covariance (ANCOVA) was conducted using oiling category as the class variable and trend as the covariable. Sampling periods for each sampling year were pooled to generate an annual trend rate (Trend). Because a curvature in mean length data for the oiling categories was noted, an annual trend squared term (Trend²) was added to the model. The model was calculated in two ways: an effects model tested for differences among means and a means model provided estimates of the means. The effects model tested for differences among oiling categories. The means model provided estimates for all of the parameters (tested from zero, not from each other) and standard errors.

Various correlations were sought using Spearman rank correlation analysis. Spearman is a nonparametric measure of rank correlation between variables. L. irrorata density and Uca burrow density were analyzed for correlation with TPH, total above ground biomass, S. alterniflora aboveground biomass, and 0–6 cm belowground biomass.

Density data were standardized to individuals m⁻² for presentation. Statistical significance was defined as p ≤ 0.05 and p-values are reported to 2 decimal places. In cases where p = 0.00, the p-value is reported as p < 0.01. All error terms are expressed as standard error (SE) with n = 7 for each oiling category.

TPH

Sites were well differentiated by TPH supporting the establishment of the three oiling categories (Table 1). Mean TPH through all sampling periods was 0.3 ± 0.1 at RF, 6.0 ± 1.9 at MD, and 106.3 ± 50.1 mg g⁻¹ at HV sites. TPH at the HV sites declined from approximately 500 mg g⁻¹ at 9 mo after the spill as reported in Lin et al. (2016), to about 50–150 mg g⁻¹ 66 mo after the spill. Similarly, MD sites declined from approximately 70 mg g⁻¹ 9 mo after the spill to 2–10 mg g⁻¹ 66 mo after the spill.

Vegetation

Lin et al. (2016) detail vegetation recovery to 42 mo after the spill; Figure 2 extends their data for total aboveground biomass to 66 mo after the spill. The total aboveground biomass did not recover at HV relative to RF and MD sites through 66 mo after the spill. From 30 to 66 mo after the spill, S. alterniflora comprised 53.28% ± 4.19% of the total aboveground biomass at RF; 66.39% ± 3.95%, at MD; and 100%, at HV sites. Lin et al. (2016) found that S. alterniflora aboveground biomass recovered at HV sites from 24 to 36 mo after the spill, and that aboveground biomass of S. alterniflora did not differ at MD across time when compared to RF sites. Aboveground biomass of J. roemerianus recovered at MD sites within 30 mo after the spill and remained at levels similar to RF sites 66 mo after the spill. Aboveground biomass of J. roemerianus at HV sites did not recover 66 mo after the spill (Fig. 2). Lin et al. (2016) also provide an analysis of live belowground biomass (g m⁻²). Live belowground biomass in the uppermost 6 cm was similar at RF and MD sites after 24 mo. Belowground biomass in this interval remained lower at HV compared to the MD and RF sites through 66 mo.

Uca species, burrow density, and burrow size

Uca density varied from 5 to 30 burrows m⁻² (Fig. 3) and was significantly influenced by sampling period (p < 0.01), but without influence by oiling intensity (p = 0.46), although the interaction between sampling period and oiling intensity was observed (p = 0.16). Tukey’s test did not detect significant differences among sampling periods. However, the potential for an ongoing impact is indicated by large differences between both HV and MD sites compared to RF sites at 30 and 48 mo after the spill (Fig. 3). Uca burrow diameter (Fig. 4) varied from <5 to 30 mm. Diameter significantly varied by sampling period (p < 0.01), but not oiling intensity (p = 0.90) and without interaction (p = 0.91). Tukey’s test again indicated no significant differences among sampling period, however a potential ongoing impact is indicated by a large difference in burrow diameter between both HV and MD compared to RF sites at 30 and 54 mo after the spill.

Uca longisignalis was the most common fiddler crab species overall and was observed at all sites and sampling times. However, we observed Uca spinicarpa only at HV sites and only early in our sampling period, 30 and 36 mo post-spill. The marsh platform at the sites where U. spinicarpa was observed was largely unvegetated with a residual oil crust on the surface.

Littoraria irrorata density

L. irrorata total density (all length classes combined, Fig. 5) was influenced by sampling period (p < 0.01), but not oiling intensity (p = 0.08) and without interaction (p = 0.49). Similarly, the densities of adults, subadults, and juveniles of L. irrorata all varied across all sampling periods, each with p < 0.01, but not by oiling intensity (p = 0.08, 0.08, and 0.31, respectively) and without interaction (p = 0.38, 0.20, and 0.06, respectively). There was a noticeable trend of increasing total population density from 30 to 66 mo after the spill, except for a decrease at the HV sites at 36 mo after the spill. Total population density at MD increased over time at a faster rate than either RF or HV sites, which increased at similar rates. Total population density at MD remained higher than RF or HV sites even at 66 mo after the spill. Beginning at 48 mo after the spill, densities at all sampling sites (most noticeably at MD sites) generally appeared to vary seasonally with higher density in the spring and lower density in the fall. Although not significant at the p ≤ 0.05 level, densities at MD sites strongly trended higher compared to the RF and HV sites at 40 and 48 mo after the spill (Fig. 5).

Littoraria irrorata shell length

Mean L. irrorata shell length (Fig. 6) varied from ∼16–21 at RF sites, ∼16–19 at MD sites, and ∼12–18 mm at HV sites. Mean shell length was influenced by sampling period (p < 0.01), oiling intensity (p = 0.05), and the interaction of oiling and sampling period (p < 0.01). A general trend of increasing mean shell size over time for all of the oiling levels is noticeable when data from 30 to 40 mo are compared with data from 48 to 66 mo. This observation was most distinct at RF sites.

Trends over time were further investigated using ANCOVA with oiling category as the class variable and “Trend” as the covariable. Results reveal significant differences among oiling categories and in linear and quadratic trends for each oiling category (Table 2) providing strong support that trends over time were significant.

A separate means model was also calculated to determine if these trends differed among oiling categories, and results for intercepts and slopes were also significant, with the exception of RF sites (Table 3). Average shell length at RF sites gradually increased over the sampling period (Fig. 7). These trend lines indicate that snails at MD may have attained a similar mean shell length as RF sites between 42 and 48 mo after the spill; thereafter, however, mean shell length at MD sites decreased to the last sampling event while the mean shell length at RF sites increased. Mean shell length at HV sites increased over the study to a maximum shell length between 48 and 54 mo after the spill, but also decreased from 54 to 66 mo.

To better evaluate changes in mean shell length, the relative frequencies of L. irrorata juvenile (<6 mm), subadult (6–13 mm), smaller adult (13–20 mm), and larger adult (>20 mm) age classes were calculated by pooling size classes into yearly intervals from 30 to 66 mo after the spill (Fig. 8). The relative frequency of subadults and smaller and larger adults at RF and MD sites remained nearly equally high through 36 mo after the spill. A large proportion of smaller adults occurred at all of the oiling categories from mo 36 to 48 after the spill. However, the relative frequency of larger adults at MD and HV sites decreased after mo 48. The differences in the proportion of larger adults became more apparent from mo 48 to 60 after the spill when RF sites had a higher relative frequency than MD sites, and MD had a higher relative frequency than HV sites. The proportion of larger adults at HV sites remained nearly the same, less than 15% throughout the study, whereas larger adults comprised approximately 50% of the population at RF sites.

Correlations

Neither L. irrorata nor Uca burrow density correlated with TPH (p = 0.46 and 0.58, respectively). L. irrorata density correlated with total aboveground biomass (p = 0.01), S. alterniflora aboveground biomass (p < 0.01), and 0–6 cm belowground biomass (p < 0.01). Uca burrow density and total aboveground biomass (p = 0.74), S. alterniflora aboveground biomass (p = 0.12), and total 0–6 cm below ground biomass (p = 0.70) were not correlated.

Uca

Zengel et al. (2016b) synthesized data from five post-DWH oil spill studies by comparing oiled sites to reference sites in each study. This analysis included our data through 2014 (54 mo after the oil spill). Species composition, burrow density, and burrow size were evaluated in relation to oiling. Zengel et al. (2016b) concluded that all three metrics indicated that the DWH event negatively affected fiddler crab populations in Barataria Bay; crabs at oiled sites were less abundant through 48 mo after the spill, recovered in size within 24 mo after the spill, and differed in species composition from crabs at reference sites. Although our data were included, our results differ in some subtle ways from other studies. We found that species composition differed among RF and HV sites at 30 mo (2.5 years) after the spill, but that recovery of species composition occurred by about ∼36 mo (three years). However, our results agree with the conclusions of Zengel et al. (2014) that the loss of vegetation and change in sediments led to the reduced dominance by U. longisignalis and an increased presence of U. spinicarpa at the HV sites. Similar to the meta-analysis, we found that crab size recovered before 30 mo after the spill. The estimate from our sites for recovery of burrow density indicates a more rapid recovery than the overall estimates based on Zengel’s et al. (2016b) analysis which included sites with heavier oiling than observed in our study.

Littoraria irrorata

Zengel et al. (2016c) provides a synthesis of oiling effects on L. irrorata through mo 60 after the spill. This meta-analysis concluded that overall density at heavily oiled sites was lower than reference sites through 60 mo after the spill, although the difference decreased between mo 12 and 36 without achieving recovery. Our data reflect a slightly lower density at HV sites through mo 48 and recovery by mo 60. Zengel et al. (2016c) showed that mean L. irrorata shell length tended to be shorter at oiled sites than at reference sites through mo 12, recovery by mo 24, and a decline relative to reference sites from mo 24 to 48. We observed a similar decline after mo 36. Our data also agree with Zengel et al. (2016c) in that a greater relative proportion of smaller adults and sub-adults, and fewer large adults, were found at oiled compared to reference sites. Zengel et al. (2016c) also found a lower relative proportion of juveniles at oiled compared to reference sites, indicating low recruitment or poor survival at oiled sites. This meta-analysis also demonstrates a depression in the mean shell length in the reference L. irrorata population in Barataria Bay in 2012 (30 mo) and 2013 (42 mo), as observed in our study, and recovery by 2015 (66 mo) (Fig. 7). This depression throughout Barataria Bay could be the result of the oil remaining on the water from the spill and cleanup activities (see below). The mean shell length (Table 3 in Zengel et al., 2016c) at reference sites close to the marsh edge of all studies ranged from 14.2 to 19.6 mm. L. irrorata at our RF sites reached that length by mo 66. We found that recovery at shell length had not occurred at MD and HV sites 66 mo after the spill.

Zengel et al. (2016a) (also found in Baker, Steinhoff & Fricano, 2017; National Ocean and Atmospheric Administration, 2016; Powers & Scyphers, 2016) proposed a definition of recovery of L. irrorata as attaining densities and shell-size distributions similar to reference conditions, and predicted that recovery would take at least 36–60 mo after the oiling, and perhaps as long as 10 years or more for the largest, oldest snails, once oiling and habitat conditions were suitable, to support normal levels of periwinkle recruitment, immigration, survival, and growth in impacted areas. Although initial increases in shell size were observed at MD and HV sites, we found a decline in the relative abundance of larger adults occurred in mo 36–48 and 48–60 after the spill. If the rate of recovery was continuous, one would expect that these larger adults should attain a high relative abundance either by growth or immigration by 60 mo after the spill. This discrepancy could mean that either L. irrorata individuals are not surviving to these larger sizes, are growing more slowly than expected, or that adults leave the oiled area prior to growth to these larger sizes. Larger L. irrorata may be emigrating from the oiled areas. Vaugh & Fisher (1992) indicated a limited overall range of movement for L. irrorata; however, Silliman et al. (2005) indicate that movement of 2 m or more could occur when motivated by lack of habitat and food source or potential predation. TPH levels in the soil at HV sites remained elevated over time (Table 1). These levels have the potential to cause chronic toxicity effects on L. irrorata because snails graze both on the marsh surface and on S. alterniflora plants. Certainly, elevated TPH could indicate the presence of polycyclic aromatic hydrocarbons (PAHs) on the soil surface (Turner et al., 2014) where L. irrorata feeds upon the microalgae and other organisms. Mohammed (2015) studied PAH in the tissue of S. alterniflora and L. irrorata in lightly and heavily oiled sites and hypothesized that periwinkle grazing on S. alterniflora, where PAHs were incorporated into the plant tissues, could present an additional exposure pathway for the snails. Oiling effects on habitat quality or food resources that are slow to recover but specifically impact adult periwinkles may also explain the decline in larger adults at oiled sites.

The extent of shoreline oiling in northern Barataria Bay (Michel et al., 2013; Nixon et al., 2016) and the extent and duration of oil in and on the water estimated by the Environmental Response Management Application (ERMA) (https://response.restoration.noaa.gov/maps-and-spatial-data/environmental-response-management-application-erm; MacDonald et al., 2015) indicates that populations of L. irrorata were likely affected in their ability to carry out their life cycle on the marsh and in the water column (as larvae) at all sites (including our reference sites) for some period following the spill (Pennings et al., 2016). Even at 30 mo after the spill when this study began, the mean shell size of L. irrorata was smaller than would be predicted by regional averages at all oiling levels, including reference sites, and the recovery process was ongoing at all sites, as indicated by our analysis (Fig. 7) (also see Fleeger et al., 2017 for a similar discussion of infaunal responses).

The role of Spartina alterniflora in macrofauna recovery

We found a significant correlation of L. irrorata density with S. alterniflora aboveground biomass, and L. irrorata was found to have achieved and maintained higher density of individuals more rapidly over the study period at MD compared to RF and HV sites. S. alterniflora aboveground biomass and L. irrorata density both were equivalent at MD sites compared to RF sites 24–36 years after the spill; however, HV sites did not recover until about 60 mo or more after the spill (Fig. 5). Many studies (Kiehn & Morris, 2009; Silliman & Zieman, 2001; Silliman & Bertness, 2002; Silliman & Newell, 2003; Silliman et al., 2005) have noted that periwinkle density is positively correlated with S. alterniflora aboveground biomass. Stagg & Mendelssohn (2012) found that L. irrorata growth, survival, and productivity were positively correlated to increasing S. alterniflora canopy cover in restored marshes. Lin et al. (2016) showed that S. alterniflora achieved and maintained higher stem density at MD sites compared to RF and HV sites after 40 mo (Fig. 6B, Lin et al., 2016), further supporting a strong relationship between S. alterniflora and L. irrorata, which also achieved higher density at MD sites after 40 mo (Fig. 5).

Although we found no correlation across our sampling period, the observed pattern of recovery of Uca burrow density and species composition also tracked recovery time of live aboveground biomass of S. alterniflora (24–36 mo) noted at HV sites by Lin et al. (2016). Lin et al. (2016) also found that recovery of J. roemerianus was much slower than that of S. alterniflora at HV sites, suggesting that S. alterniflora was critical to recovery processes for both fiddler crabs and periwinkles. Lin et al. (2016) noted that MD sites were equivalent in live aboveground biomass of both S. alterniflora and J. roemerianus and live belowground biomass by 30 mo after the spill; whereas, the HV sites had no aboveground biomass of J. roemerianus and reduced recovery of live belowground biomass at 42 mo–60 mo after the oil spill. The 24–36 mo recovery of fiddler crab burrow density also parallels recovery of the most abundant meiofauna groups at the HV sites (Fleeger et al., 2015). Live belowground biomass of roots and rhizomes has proven to be very slow to recover at the HV sites, although this measure of marsh health did not appear to influence fiddler crab density or size or periwinkle density. Thus, joint, long-term studies of plants and benthos (this study and others) provides a better picture of recovery than individual ecosystem components. In conclusion, we found that the marsh as a whole did not recover within 66 mo after oiling (including periwinkle population structure in this study, and some vegetation and infaunal components; Lin et al., 2016; Fleeger et al., 2017). Our team will continue to follow recovery with proposed additional study with special interest in examining the relationships of recovery among macroinvertebrates and other marsh ecosystem components, including the ongoing recovery of vegetation and other biota.