Drought-Induced, Punctuated Loss of Freshwater Mussels Alters Ecosystem Function Across Temporal Scales

DuBose, Traci P.; Atkinson, Carla L.; Vaughn, Caryn C.; Golladay, Stephen W.

doi:10.3389/fevo.2019.00274

ORIGINAL RESEARCH article

Front. Ecol. Evol., 18 July 2019

Sec. Population, Community, and Ecosystem Dynamics

Volume 7 - 2019 | https://doi.org/10.3389/fevo.2019.00274

This article is part of the Research Topic Death and Decomposition in Aquatic Ecosystems View all 11 articles

Drought-Induced, Punctuated Loss of Freshwater Mussels Alters Ecosystem Function Across Temporal Scales

$\nTraci P. DuBose$ Traci P. DuBose¹^*

Carla L. Atkinson²

Caryn C. Vaughn¹

Stephen W. Golladay³

¹Oklahoma Biological Survey & Department of Biology, University of Oklahoma, Norman, OK, United States
²Department of Biological Sciences, University of Alabama, Tuscaloosa, AL, United States
³Jones Center at Ichauway, Newton, GA, United States

Punctuated, mass mortality events are increasing for many animal taxa and are often related to climatic extremes such as drought. Freshwater mussels are experiencing increased mass mortality events linked to hydrologic drought. Because mussels play important functional roles in rivers it is important to understand the ecosystem effects of these die-offs. Here, we address how mass mortality events of mussels caused by drought may impact stream ecosystem function. We first present a conceptual model, based on the literature, of how mussel mass mortality should affect different ecosystem functions across various ecological time scales, from hours to decades. Next, we highlight two case studies of drought-linked, mussel-mass mortality events from rivers in the southern U.S. We then present the results of an experiment we performed quantifying the ecosystem effects of a punctuated mussel die-off. Finally, we combine our experimental results with field data from a recent mussel die-off to predict how mussel losses will influence ecosystem function. Based on the presented case studies, our mesocosm experiment, and our extrapolated nutrient pulse due to a mussel die-off, we conclude that stream ecosystems are extensively altered following mussel mass mortality events. Mussel loss is governed by drought severity, location within the river network, and species-specific drought tolerances. In the short term, decomposing carrion from mussel die-offs releases a large pulse of nutrients into the water which stimulates food web productivity. In the long term, the overall loss of mussel biomass, and the loss of functional traits as more sensitive species decline, leads to decreases in ecosystem function which may take decades to recover. Drought and human demand for water will make mussel die-offs more likely in the future and it is unlikely that drought sensitive species will recover without changes in water management and restoration of populations through mussel propagation. Our research provides an example of how the loss of an abundant, long-lived organism has cascading, and long-term impacts on ecosystems.

Introduction

Resource pulses are episodes of increased resource availability in space and time that are relatively rare, of large magnitude, and usually of short duration (Yang et al., 2008). These pulses are widespread and often result from climatic and environmental factors. Resource pulses can result from the mass die-offs of animals, such as 17-year cicadas, spawning salmon, and even wildebeest, and are increasingly recognized as important components of ecosystem function (Yang et al., 2008; Subalusky et al., 2017). Mass mortality events, or die-offs, are increasing in frequency across most taxa (Fey et al., 2015), thus it is important that we understand how these events affect ecosystem function (Baruzzi et al., 2018). In freshwater systems, unionid mussels play important structural and functional roles and are also experiencing increasing mass mortality events globally (Lydeard et al., 2004; Wenger et al., 2018) that are often linked to climatic events such as drought (Vaughn et al., 2015). Thus, they are a useful system for investigating the ecosystem effects of mass mortality events, particularly as related to environmental change (Fey et al., 2015).

Here, we address how mass mortality events of unionid mussels caused by drought may impact stream ecosystem function short-term and long-term. We first present a conceptual model, based on the literature, of how mussel mass mortality should affect different ecosystem functions across various ecological time scales, from hours to decades. Next, we highlight two case studies of drought-linked, mussel-mass mortality events from rivers in the southern U.S. We then present the results of an experiment we performed quantifying the ecosystem effects of a punctuated mussel die-off. Finally, we combine our experimental results with field data from a recent mussel die-off to predict how mussel losses will influence ecosystem function.

How do Mussel Mass Mortality Events Impact Ecosystem Function? A Conceptual Model

Freshwater mussels (order Unionida, hereafter mussels or unionids) are sedentary mollusks that live burrowed in stream sediments where they filter the water and transfer energy and nutrients from the water column to benthos. Nutrients excreted and biodeposited by mussels stimulate instream microbial, primary, and secondary production and are even exported to riparian areas (Allen et al., 2012; Vaughn, 2018). Mussels provide habitat for other organisms through the biogenic structure of their shells and by changing hydrodynamic conditions at the sediment-water interface (Sansom et al., 2018b). They are long-lived (6–100 years) with high native biodiversity in eastern North America (Williams et al., 1993), often live in high density, multi-species aggregations (hereafter mussel beds) that can persist in rivers for many decades and can make up most of the invertebrate biomass in many perennial rivers (Sansom et al., 2018a). Mass mortality of mussels has been linked to increasing drought, either from emersion (Atkinson et al., 2014) or from low dissolved oxygen and high temperatures associated with decreased water volume (Gagnon et al., 2004), as drying in streams is often accompanied by increased water temperatures and diel oxygen shifts (Mosley, 2015). How mussels respond to drought conditions depends on individual species' physiological tolerances, drought severity, and abiotic conditions (Gagnon et al., 2004; Golladay et al., 2004; Haag and Warren, 2008; Gough et al., 2012). When mass mortality of mussels occurs, their loss should influence stream ecosystem functions in a variety of ways across ecological time scales (Figure 1). We highlight these predicted effects below.

FIGURE 1

Figure 1. Conceptual model depicting ecosystem function shift due to unionids at short and long time scales. Pre mass mortality information based on 25°C, downstream values for 1991 in Vaughn et al. (2015) for the biofiltration and community structure data and on Atkinson and Vaughn (2015) for downstream nutrient excretion and storage information. *Loss of nutrient capacitance dependent on mussel recruitment. Time increases to the left and the impact on mussels of four main ecosystem functions is described and depicted. Capital letters refer to distinct processes in the model, while lower case letters refer to sources: (a) Vaughn et al. (2015), (b) Atkinson and Vaughn (2015), (c) McDowell et al. (2017), (d) Gagnon et al. (2004), (e) Cherry et al. (2005), (f) Cherry et al. (2005), (g) Atkinson et al. (2014), (h) Atkinson et al. (2018), (i) Ilarri et al. (2018).

Mussel Mortality Feedback Loop

Mussel species vary in their physiological tolerance and response to stress (Spooner and Vaughn, 2008). Species that are sensitive to low oxygen (hypoxia sensitive) or higher temperatures (thermally sensitive) are less likely to survive during a hydrologic drought, thus thermally tolerant and/or hypoxia tolerant mussels become the dominant species within the assemblage (Gagnon et al., 2004; Atkinson et al., 2014). Mussel soft tissue can decay within seven days, as shown below. Decomposing soft tissue releases a pulse of nitrogen and phosphorus into the water column and interstitial spaces (Figure 1A; Cherry et al., 2005; Atkinson et al., 2014). Depending on stream discharge, this nutrient pulse moves downstream over a few hours/days. Pore water can retain high nutrient concentrations longer, potentially exposing burrowing unionids to lethal nutrient concentrations (Cooper et al., 2005; Gough et al., 2012). After the loss of mussels, shifts in algal production, and turbidity are driven by stream discharge. If the stream becomes intermittent (a series of drying pools), turbidity will likely decrease due to increased sedimentation while algal blooms will likely form in stagnant areas (Mosley, 2015). During intermittence, the combination of reduced biofiltration and the release of nutrients into the water from mussel soft tissue decay encourages large algal blooms (Gagnon et al., 2004), which leads to high respiration rates at night, further reducing dissolved oxygen concentrations and stressing the remaining mussels, leading to additional mortality (Figure 1B). This cycle can also exacerbate feedback among deaths within mussel beds; algal blooms cause mortality in remaining hypoxia-intolerant bivalves, worsening algal blooms, and depressing dissolved oxygen, further stressing and eventually killing hypoxia-tolerant mussels. If the stream remains perennial, turbidity will increase through the addition of suspended solids from upstream and algal blooms become less likely. Algal blooms and/or increased turbidity can persist because of a reduction in biofiltration by freshwater mussels (Figure 1C).

Reduced Biofiltration

Immediately after a mass mortality event, biofiltration is greatly reduced in part to both residual stress on the remaining living mussels and biomass loss from mussel mortality (Vaughn et al., 2015). Biofiltration by the remaining mussels will gradually increase within the following days, but is likely to remain low until mussel biomass is replaced, unlikely for at least a decade (Figure 1E). Reduced biofiltration drastically increases the time required for the remaining mussels to filter a given amount of water, reducing material exchange between the water column and benthos (Baustian et al., 2014; Vaughn et al., 2015).

Reduced Nutrient Capacitance and Storage

Short-term nutrient storage in mussel soft tissue is greatly reduced through decomposition (Atkinson et al., 2014). As mussels filter feed, they act as “nutrient capacitors,” accumulating, storing, and releasing energy and nutrients (carbon, nitrogen, and phosphorus) at different rates based on their age and species' traits (Strayer, 2014; Atkinson et al., 2018). Following mass-mortality, remaining mussel assemblages have lower abundance, age diversity, and species diversity, reducing their ability to filter seston and excrete nutrients. This reduced nutrient capacitance may result in longer nutrient spirals and more downstream transport of nutrients, likely due to an increase in nutrient uptake length (Figure 1D; Atkinson et al., 2014). Shells of deceased mussels lose ~50% mass by 15 years through mechanical and chemical dissolution, which reduces the nutrient storage and shell habitat within the mussel bed as particle size becomes more homogeneous (Figure 1E; Atkinson et al., 2018).

Changes in Habitat Provided by Shells

Live mussels and their spent shells physically modify the environment in streams, providing unique habitat for other organisms. Tissue decay potentially creates interstitial spaces within the substrate, which can be used by both macroinvertebrates and fish. Shells vary in shape and size across species and age and can accumulate in the sediment at different rates; thus, shell habitat can harbor variable macroinvertebrate communities depending on the shells species of origin (Bódis et al., 2014). While shells represent hard surfaces for macroinvertebrates, they dissolve over time. Shell dissolution is fastest in flowing waters with low calcium carbonate concentrations and thin, small shells dissolve faster than thicker, larger shells (Strayer and Malcom, 2007; Ilarri et al., 2015, 2019). While bivalve soft tissue decomposes quickly, shells persist for many decades (5–30 years), providing habitat for other stream biota (Strayer and Malcom, 2007; Ilarri et al., 2015; Atkinson et al., 2018). Over time, the benthos will be more homogenous as old shells dissolve and new relic shells are produced by a less diverse mussel assemblage ultimately altering benthic microhabitat characteristics and macroinvertebrate community structure (Figure 1E; Ilarri et al., 2018).

Shifts in Community Composition and Ecosystem Function

During droughts, species sensitive to low oxygen (hypoxia sensitive) or higher temperatures (thermally sensitive) face greater risk of mass mortality leading to differential survival resulting in tolerant species becoming dominant within an assemblage, changing community structure (Gagnon et al., 2004; Atkinson et al., 2014). Surviving mussels may contribute to population recovery if conditions are suitable for reproduction. Most mussels have an ectoparasitic larval phase that requires a host fish (Barnhart et al., 2008). Drought concentrates fish into drying pools as they to attempt to escape harsh conditions or die due to increased biotic and abiotic stressors (Matthews and Marsh-Matthews, 2003; Lennox et al., 2019). While mussel reproduction is unlikely limited by host density (Haag and Stoeckel, 2015), different mussel species exhibit different host specificity and infection phenology (Barnhart et al., 2008). Thus, predicting the recruitment success is difficult due to unionid's unique life histories. As it takes mussels anywhere from 9 months to 10 years to reach sexual maturity, with most mussels reaching maturity around 4 years old (Haag, 2012), in die-off affected areas, biomass is unlikely to rebound for at least a decade. Future mussel assemblage structure is dependent on the surviving mussel assemblage, the surviving fish assemblage, and the recurrence frequency of droughts. If drought frequency decreases or remains constant, the mussel community could return to its former, pre-drought structure if no mussel species were extirpated from the river basin. If droughts increase in frequency and severity as projected in many regions (Palmer et al., 2008; USGCRP., 2017), we anticipate the mussel community will become dominated by tolerant mussel species (Figure 1E).

Mussel species with different temperature tolerances have different, temperature-dependent biofiltration and nutrient excretion rates. Thus, when the proportion of thermally sensitive vs. tolerant species in a mussel assemblage changes, this can impact ecosystem function (Spooner and Vaughn, 2008; Vaughn et al., 2015; Atkinson et al., 2018). For example, in rivers in southern Oklahoma, Actinonaias ligamentina is a thermally sensitive species with higher filtration and nutrient excretion rates at summer temperatures than other mussels in the assemblage. Because of its temperature intolerance, it also has a higher mortality rate during drought than other species. Thus, when this species is lost, the overall biofiltration (Figure 1C) and nutrient recycling capacity (Figure 1F) of the community is decreased for an extended time period, even if the biomass of other species remains stable (Vaughn et al., 2015). Further, mussel species excrete at different N:P ratios, and losses can also lead to shifts in excretion stoichiometry (Figure 1F; Atkinson et al., 2018). These changes can cascade through the food web impacting algal, macroinvertebrate (Novais et al., 2017), and even fish dynamics.

Two Case Studies of Drought-Driven Mussel Losses in the Southern United States

Drought-induced mussel mass mortality events have been documented for two diverse, well-studied river systems in the southern U.S., the Lower Flint River in Georgia and the Kiamichi-Little River system in Oklahoma (Figure 2). These case-studies allow for a deeper understanding of how drought affects mussel assemblages, the subsequent changes in stream ecosystems, and the potential recovery time for mussel assemblages.

FIGURE 2

Figure 2. Map of case study locations within the southern continental USA. In the bottom panel, large streams (stream order ≥8) are depicted in light blue to provide context while the focal river rivers are found within the red boxes. The letters next to the river basins indicate which panel corresponds to those river basins. Southeastern Oklahoma (A) contains three rivers: the Kiamichi, Little, and the Mountain Fork. K2 and K3 are sites on the Kiamichi discussed in the ‘short and long-term nutrient releases following bivalve mortality' section of this manuscript. The lower Flint River (B) is in southwestern Georgia; rivers discussed in the case study section are named. Code from the R package hydroMap (DeCicoo and Blodgett, 2017).

Lower Flint River, Georgia

The lower Flint River of southwestern Georgia experienced an extended period of below normal rainfall from 1999 until 2013. This period included three multiyear droughts that were classified as severe/exceptional. The summer of 2000 was particularly disastrous for unionids as streams in the region experienced unprecedented low flows and transitioned from perennial to intermittent. Forty-six historically-species-rich sites in lower Flint tributaries were surveyed in 1999 prior to drought onset, and a subset of these were resampled in 2000 (Gagnon et al., 2004). Stream drying had not been previously observed but became common during subsequent growing seasons (Rugel et al., 2012). Stream flow was essential for maintaining dissolved oxygen concentrations within the tolerances of freshwater mussels; dissolved oxygen between 5 and 3 mg L⁻¹ resulted in 24% mortality and when dissolved oxygen fell below 3 mg L⁻¹ up to 76% of mussels died (Gagnon et al., 2004; Golladay et al., 2004). These degraded physicochemical conditions differentially impacted species within the mussel communities. Riffle specialists in medium-sized streams suffered the highest mortality, while drought-tolerant, small stream species and species in larger tributaries whose habitat was buffered from drought conditions fared better (Gagnon et al., 2004). In medium-sized streams, community composition shifted toward more generalist species from riffle specialists (Gagnon et al., 2004).

The lower Flint River has been surveyed since this drought. Through 2013, there was little evidence of recovery from mortality associated with the 1999–2001 drought (Smith et al., 2015). During most summers since the initial drought, conditions were stressful, likely preventing the reproduction of surviving unionids. Reproduction of mussel populations, as evident from observation of juvenile mussels, was not apparent until rainfall approached average levels (2013–2015) (Smith et al., 2015). The extended period of below normal rainfall (1999–2013) and subsequent below normal stream flows were likely exacerbated by anthropogenic water withdrawal; the mid reaches of tributaries of the lower Flint cross the Dougherty Plain physiographic district, which is a recharge area for a heavily developed agricultural water source (Golladay et al., 2004). This case study provides evidence that mussel biomass might not recover from mass mortality events for over a decade and that anthropogenic and climate alterations can alter stream benthic communities for extended periods.

Kiamichi and Little Rivers, Oklahoma

The Kiamichi and Little Rivers in southeastern Oklahoma are adjacent, major tributaries to the Red River. This region experienced a period of exceptional drought during 2011–12 where the Kiamichi River experienced 84 days of no flow (defined as discharge < 0.01 m³ s⁻¹) and 36 weeks of extreme low flow, defined as flows below the 10th percentile of flow frequency (Atkinson et al., 2014). The Little River, and its major tributary the Mountain Fork River, experienced 39 and 40 weeks of extreme low flow, respectively (Atkinson et al., 2014).

These severe drought conditions led to a mass mortality event as mussels became isolated in shallow drying pools or emersed. Mussel losses and their effects on ecosystem function were documented in two related studies. Atkinson et al. (2014) sampled mussels at three sites before a drought (2010) and at the end of the drought (2012) and assessed changes in mussel abundance and mussel-provided nutrient cycling and storage. Sixty percent of unionids died during the drought, but thermally sensitive species had a higher mortality rate, resulting in a community shift toward more thermally tolerant species. In the second study, Vaughn et al. (2015) compared mussel biomass and ecosystem services (biofiltration, nutrient cycling, and nutrient storage) at four sites in the Kiamichi River across several decades (1991, 2004, 2011). 1991 was a wet period and 2004 and 2011 were drought periods with significant mussel losses. They found that mussel biomass decreased over 60% across these sites and that ecosystem function losses mirrored the biomass losses. The sites experienced mussel biomass losses of ~28% and corresponding declines in nitrogen recycling (22%), phosphorus recycling (15%), and ~30% declines in areal storage of nitrogen and phosphorus (Vaughn et al., 2015).

An Experiment Quantifying the Effects of a Mussel Die-off on Ecosystem Function

While case studies have documented how native mussel communities change after mass mortality events, no studies have experimentally demonstrated how these losses impact ecosystem function. With our conceptual model and case studies in mind, we designed an experiment to measure ecosystem function changes that occur following a punctuated mussel mass mortality event. We conducted a mesocosm experiment at the University of Oklahoma Biological Station in the summer of 2018 where we induced mussel mortality and measured effects on ecosystem structure (water column nutrient concentrations and algal abundance) and ecosystem function (decomposition rates and ecosystem metabolism) over time. We predicted that decaying mussel tissue would increase nutrient concentrations, which would stimulate both algal growth and microbial respiration and decomposition.

Mesocosm Experiment Methods

We used 18, 1.52 m diameter, 946 L circular tanks to simulate drying stream pools. Each mesocosm (tank) was lined with ~15 cm of gravel (1:1 ratio of 10 and 38 mm diameter gravel). We had 9 control mesocosms with no mussels and 9 mesocosms containing 31 mussel individuals, to replicate a natural mussel community in the region. Each mussel treatment mesocosm contained 13 Actinonaias ligamentina, 9 Cyclonaias pustulosa, 5 Amblema plicata, two Tritogonia verrucosa, one Lampsilis cardium, and one Plectomerus dombeyanus. This represented a low, but natural density of mussels (11.9 mussels/m²) and reflected the freshwater mussel community of the upper Kiamichi River (Atkinson et al., 2012).

We describe sampling events and present our results relative to the day of the mussel mass mortality event: negative values indicate days before and positive values days after the mussels died. We filled mesocosms with water 12 days (day−45; see Supplementary Table 1 for sampling dates and measurements) before adding mussels to allow the mesocosms to be naturally colonized by algae and macroinvertebrates. On day−33, we added mussels. On day−14, we introduced 10 largemouth bass (Micropterus salmoides; mean standard length = 95 mm, SD = 16 mm) to simulate how fish are concentrated in drying pools during early periods of drought. We removed the fish on day−2 to simulate their movement downstream during the drought as drying pools became too stressful for them (Magoulick and Kobza, 2003). Fish may have impacted water column nutrients on day−3 but had little impact on mesocosm nutrient concentrations and algal abundance (see Results below). We induced a punctuated mass mortality event on 2 July 2018 (day 0) by sacrificing the mussels in 5 of the 9 mussel mesocosms, while maintaining the 9 non-mussel controls. To produce the carrion for this stage of this experiment, we cut the adductor muscles of 155 mussels. We returned the mussel carrion to the mesocosms to allow for natural decomposition of soft tissue.

We sampled mesocosms 3 times before and 4 times after the mass mortality event, resulting in the following sampling days:−20,−15,−3, 4, 11, 25, and 39. On each sampling day we measured dissolved oxygen (DO), conductivity, and water temperature at midday (Supplementary Table 2). To measure water column phosphorus and nitrogen concentrations, we collected 20 mL filtered water samples (Grade F, 0.7 μm pore size, Sterlitech Kent, Washington) and froze them for subsequent nutrient analysis. We lost nutrient samples from day 11, thus we resampled the mesocosms on day 18. We filtered two water samples (Grade A, 1.6 μm pore size, Sterlitech Kent, Washington) from each mesocosm to quantify water column chlorophyll a. Filters were frozen for later chlorophyll estimation. On day−33, we placed six 7.6 cm² clay tiles with a 27.5 mm² fritted glass disc attached with silicone (LECO cover crucible AL P 1000; GE Silicone 1^* All Purpose) on the substrate surface to allow algal colonization for sampling of benthic algal production. We removed two glass fritted discs on each sampling day and froze them for later estimation of benthic algal biomass.

We quantified soluble reactive phosphorus (SRP) with the colorimetric method (Murphy and Riley, 1962; Stainton et al., 1974; EPA Method 365.3) and ammonium (NH₄-N) using the phenol method (5.2.6 EPA Method 350.1; ASTM., 2012) for the filtered water samples. To measure chlorophyll a concentration, we cold-extracted water column (filters) and benthic (fritted discs) samples with acetone and measured the extractant spectrophotometrically with a correction for pheophyton (ASTM., 2012).

We followed Tank et al. (2017) to measure ecosystem metabolism as gross oxygen production on days 4, 11, 25, and 39. We quantified ecosystem metabolism by measuring dissolved oxygen production and respiration in light and dark cycles, respectively, on the glass fritted discs in 50 mL centrifuge tubes. We measured dissolved oxygen (HACH HQ40d multiple parameter meter, Loveland, Colorado) to estimate initial oxygen concentrations. We placed fritted disks in centrifuge tubes filled with the respective mesocosm's water and sealed the tubes. After allowing the discs to metabolize for an average of 1.75 h (SD = 0.23 h) in a common mesocosm, we re-measured the dissolved oxygen within the tubes. We then removed the water and repeated the filling process with the same tube/glass fritted disc pair. After filling the tubes, we immediately placed them in the dark in a common mesocosm for an average of 2.62 h (SD = 0.47 h). We then measured final dissolved oxygen and collected and froze the discs. Gross primary production (mg DO cm⁻¹ h⁻¹) was calculated from the addition of net primary productivity (difference in DO in the light cycle) and the absolute value of the respiration (DO difference in the dark cycle).

To determine the decay rates of mussel tissue, we placed the combined shell and soft tissue of each of 5 A. ligamentina in fine mesh bags (pantyhose) in each mortality treatment (original weight mean = 297.8 g, range = 77.5–461.5 g). We chose Actinonaias ligamentina because it was the most abundant mussel species within each mesocosm, is thermally sensitive, and most likely to be lost during a drought (Atkinson et al., 2014). We weighed the bags every 12 h for 4.5 days, and then daily for 10 days until the shells were empty and the weight was stable. We calculated total tissue (including both the soft tissue and the shell) decay rates and soft tissue decay rates following Strayer and Malcom (2007). Soft tissue decay rates were determined by assuming the minimum weight measurement consisted of only shell material and subtracting that measurement from each weight measurement. To examine organic matter decomposition rates, we incubated three cotton strips (8 × 25 mm) in the bottom of each mesocosm beginning on day 7 (Tiegs et al., 2013). We removed strips on days 18, 28, and 38; these dates mirror the decomposition study by Novais et al. (2017) and reflect incubation times of 11, 21, and 32 days. Strips were preserved with 85% ethanol and later dried at 40°C. We determined the tensile strength of each cotton strip using a tensiometer (Mark 10 MG100) torn at 2 cm/min following Tiegs et al. (2013). As such, the tensile strength reflects the remaining organic matter of the original cotton strip and is reported in pounds.

All statistical analyses were conducted with R Core software version 3.5.3 (R Core Team, 2019). We used mixed linear models to test for differences among our dependent variables based on the fixed factor treatment, the fixed continuous variable sampling date, the interaction between the fixed variables, and a random intercept accounting for mesocosm. We included the mesocosm as a random intercept to account for the repeated measures over time on each replicate; this allows each mesocosm to have a different starting value and accounts for mesocosm dependency. Each model was checked visually for normality and homogeneity of variance of its residuals (Zuur et al., 2009); we log₁₀ transformed water column ammonium, water column chlorophyll a, and benthic chlorophyll a to meet these assumptions. We used the function lmer() [from the R package lme4 (Bates et al., 2015)] to perform all mixed models as we had different sample sizes for each treatment: 9 control mesocosms, 5 mortality mesocosms, and 4 live mussel mesocosms. We used the function anova() to conduct a type III ANOVA with Satterthwaithe's method and obtain p-values for all models as implemented in the R package lmerTest (Kuznetsova et al., 2017). We then used Tukey post-hoc tests to conduct multiple comparisons if the null hypothesis was rejected for each dependent variable as implemented in the package emmeans (Lenth, 2018).

Mesocosm Experiment Results

Mussel Decay

Actinonaias ligamentina soft tissue and shell had an average instantaneous decay rate of −0.016 day⁻¹ across all mortality mesocosms. Within 7 days, most soft tissue had decayed within the bags. The average instantaneous decay rate of the soft tissue alone was −0.336 day⁻¹. We did not observe shell dissolution within the time frame of our experiment.

Nutrients

Following the mass mortality event, the mortality treatments had a large increase in ammonium. Ammonium (NH₄-N) was significantly higher in the mortality treatments compared to the control treatments (F_{2, 120} = 10.92, p < 0.001; Figure 3A). Sampling day was also significant in predicting ammonium amount in the system (F_{1, 120} = 14.34, p < 0.001). The interaction between treatment and sampling day was not statistically significant (F_{2, 120} = 1.17, p = 0.31). Overall, ammonium in the water column significantly increased by 94.4% directly after the mass mortality event, while ammonium in the control mesocosm increased by 9.6%, although this was highly variable (ranged from 84.1 to −62.5%). SRP generally increased during the experiment (F_{1, 103} = 6.91, p < 0.01; Figure 3B). While SRP was not significantly different between treatments (F_{2, 16} = 0.37, p = 0.70), the interaction between treatment and sampling day was statistically significant (F_{2, 103} = 4.67, p < 0.02). Between the die-off and the end of the experiment, SRP increased 38% in tanks that experienced the mussel die-off but decreased by 51% in control tanks.

FIGURE 3

Figure 3. Ambient concentrations of NH₄-N (μg L⁻¹; A), SRP (μg L⁻¹; B), water column chlorophyll a concentrations (μg L⁻¹; C) and benthic chlorophyll a concentrations (μg cm⁻²; D) in control, live mussel, and dead mussel treatments. The dashed line represents when the mass mortality event was induced in 5 live mussel mesocosms. Points represent the mean, while the lines indicate the standard deviation. Ammonium is statistically different within control and dead treatments based on a Type III ANOVA, Tukey HSD, p < 0.05. Benthic chlorophyll a was higher in live and mortality treatments than control treatments based on a Type III ANOVA, Tukey HSD, p < 0.05.

Primary Production and Ecosystem Metabolism

Mortality treatments had higher gross primary production than the control mesocosms (F_{2, 30} = 4.11, p < 0.03). Mortality treatments had higher gross primary production than live treatments, but there was not a statistical difference between the two groups (t₂₂ = 2.27, p = 0.08). Sampling day was significant in explaining gross primary production (F_{1, 54} = 77.81, p < 0.001). The interaction between the two terms was also not statistically significant (F_{2, 54} = 1.94, p = 0.16). There was not a difference in water column chlorophyll a concentration between treatments (F_{2, 19} = 0.54, p = 0.59), although water column chlorophyll did increase through time (F_{1, 105} = 12.51, p < 0.001; Figure 3C). The interaction between the two terms was also not statistically significant (F_{2, 105} = 0.76, p = 0.47). In late July and August, some control and live mesocosms experienced algal blooms, while mortality mesocosms had low water column chlorophyll a concentrations. Benthic chlorophyll a concentration was higher in mortality and live treatments than control treatments (F_{2, 19} = 9.10, p < 0.002; Figure 3D). Sampling day did not predict benthic chlorophyll a concentration (F_{1, 106} = 0.59, p = 0.45) and the interaction between treatment and sampling day was not significant (F_{2, 106} = 0.89, p = 0.41). Note that ecosystem metabolism and benthic chlorophyll a concentrations are significantly correlated as they were measured from the same glass fritted discs (R² = 0.28, p < 0.001). While these variables are correlated, one measures ecosystem structure (biomass) and the other ecosystem function (respiration).

Organic Matter Decomposition

Higher tensile strength corresponds to a higher percentage of the original remaining cotton strip; thus, higher tension indicates less decomposition. Tensile strength of the cotton strips decreased with time (F_{1, 36} = 4.95, p < 0.04) and was significantly different between treatments (F_{2, 18} = 5.60, p < 0.02; Figure 4). The mean tensile strength of an unincubated cotton strip is 65.6 lbs (SD = 2.0 lbs).

FIGURE 4

Figure 4. Tensile strength (lbs.) as a proxy for decomposition of organic matter within different treatments. Cotton strips have a starting mean tensile strength of 65.61 lbs (SD = 1.97 lbs). Points represent the mean, while the lines indicate the standard deviation. Decomposition was higher in the mortality treatments than the live treatments and is statistically significant based on Type III ANOVA, Tukey HSD, p < 0.02.

Below, we apply these results to an actual mussel die-off to extrapolate how mussel die-offs can impact nutrient cycling in river reaches following a mass mortality event.

Scaling up: Short and Long-Term Nutrient Releases Following a Mass Mortality Event

As described above, an extreme drought in the Kiamichi River in 2011 led to massive mussel mortality. We combined quantitative data on these mussel losses with nutrient release data from our mesocosm experiment to extrapolate how mussel losses impact short and long-term nutrient cycling and storage for several river reaches.

On July 30 and 31, 2011, we (CLA and CCV) sampled mussels at three sites (K2, K3 riffle, and K3 pool) severely impacted by the drought. The K2 and K3 pools were isolated pools that still contained water but were very shallow (<10 cm) and warm (K3 exceeded 40°C) (Vaughn et al., 2015). K3 contained a riffle that was completely dry. For K2 and K3 pools, we measured mussel abundance and composition by sampling ten, 0.25 m² quadrats following Vaughn et al. (1997); we identified and measured the length of both dead and live individuals. We returned live individuals to the mussel bed. For the K3 dry riffle we laid out 10 transects across the dry mussel bed and identified shells from 14 quadrats across each transect. We used species-specific, length-soft tissue dry mass regression equations (Hopper et al., 2018) to calculate the soft tissue mass of each mussel. We calculated areal nutrient pulses (μg L⁻¹ m⁻²) released to the river as a consequence of mussels dying as the product of nutrients released by decomposing mussel tissue to the water column in the mesocosm at sampling day 4 (μg L⁻¹ g⁻¹) by the biomass loss of mussels at a site on an areal basis (g m⁻²). Based on our extrapolation, this mussel die-off resulted in a large pulse of both nitrogen and phosphorus (Figure 5). This nutrient pulse is equivalent to the areal nitrogen excretion of a mussel assemblage for 20 h, phosphorus excretion for 195 h (Atkinson et al., 2018), and to the phosphorus release from dissolving shells for two years (Wenger et al., 2018). This represents a large pulse of phosphorus that likely stimulates primary production.

FIGURE 5

Figure 5. Extrapolated nutrient release from Kiamichi mussel beds after a drought-induced mass mortality event. The shell lengths from dead unionids are used to predict mussel soft tissue (A), which was then paired with nutrient data from the mesocosm experiment to predict ammonium release (μg m⁻²; B) and SRP release (μg m⁻²; C).

Freshwater mussel soft tissue represents short term storage of nutrients since the soft tissue decays quickly, while the shell represents long term nutrient storage. We wanted to determine the role of shell material as a potential long-term nutrient sink and site of nutrient release following mortality. While shell decay can be highly variable in freshwater ecosystems (Strayer and Malcom, 2007), we used the average shell decay rates from a previously published study on the Kiamichi and Little Rivers in Oklahoma (Atkinson et al., 2018) and the Sipsey River, Alabama (Atkinson, unpublished) to estimate spent shell biomass and nutrient release (C, N, and P) over time. Specifically, using the average decay rate (k = −0.053 year⁻¹) and the average shell size per site, we estimated shell biomass (g m⁻²) over an 80-year timeframe and the subsequent nutrient release, assuming a constant rate, from spent shells following the punctuated mortality event in the Kiamichi River in 2011. We expect this estimation to be realistic as the shell decay rate was measured over a year with shell from the Kiamichi River, thus accounting for how discharge, water chemistry, and season affects shell decay rates (Strayer and Malcom, 2007; Ilarri et al., 2019).

At K3, the mass of relic shells (2.4 kg m⁻²) exceeded shell mass of living mussels (1.6 kg m⁻²), while all mussels in the dry reach perished, resulting in 1.7 kg m⁻² of shell material exposed. While some shell material may have been exported due to terrestrial scavengers, we assumed that it remained in the stream channel and was submerged once flow resumed. Site K2 did not experience complete drying and the low flows did not result in as much mortality as K3 and resulted in 1.5 kg m⁻² of relic shell while 5.0 kg m⁻² was maintained in live mussels. Based on the modeling described above, shell material decay (Figure 6A) would result in a slow nutrient release at each of these sites (Figures 6B–D). For example, 413 g C m⁻², 4.8 g N m⁻², and 0.4 g P m⁻² remained in the shell of dead mussels at K3 in the reach that did not dry, which would then be slowly released by shell decayed (Figures 6B–D). The surviving, living mussels at K3 still maintained nutrients in their shell and continued to store nutrients and potentially grow and store additional material. When mortality is equal to the production of shells, shell mass is maintained at a steady state. But large-scale die-offs lead to a pulse in relic shells and lower production of shell material. This represents a long-term loss of mussel-driven nutrient storage and shell habitat within stream reaches (Wenger et al., 2018).

FIGURE 6

Figure 6. Shell decomposition (A) of spent mussel shells following the mass mortality event as a result of low flow and high temperature conditions at 3 sites in the Kiamichi River, OK and the resultant release of carbon (B), nitrogen (C), and phosphorus (D) from the shell material over time. Shell decomposition is modeled based on the shell decomposition rates and shell stoichiometry values from Atkinson et al. (2018) and the empty shells sampled from Southeastern Oklahoma in the drought of 2011.

Implications and Conclusions

Based on the presented case studies, our mesocosm experiment, and our extrapolated nutrient pulse due to a mussel die-off, we conclude that stream ecosystems are severely altered following mussel mass mortality events. Mussel loss is governed by drought severity, location within the river network, and species-specific drought tolerances. In the short term, decomposing carrion from mussel die-offs releases a large pulse of nutrients into the water which stimulates food web productivity. In the long term, the overall loss of mussel biomass, and the loss of functional traits as more sensitive species decline, leads to decreases in ecosystem function which may take decades to recover (Figure 1).

While we have frequently observed algal blooms in the field following mussel die-offs, we did not observe algal blooms within our mesocosm experiment. In our small mesocosms, the decomposition of mussel tissue likely altered the microbial community to favor heterotrophs, which potentially out-competed water column algae for available nutrients. Our extrapolation from the observed mussel die-off in the Kiamichi River predicted a large phosphorus pulse. After inducing a die-off of the invasive bivalve Corbicula, McDowell et al. (2017) observed a smaller increase in ambient phosphorus concentration than expected. They posited that algal uptake of SRP accounted for difference between the predicted increase in water column SRP and what was measured. Further exploration of the interacting factors driving algal bloom formation after mussel die-offs during drought is warranted.

The frequency and severity of hydrologic drought is predicted to increase in the southcentral and southeastern U.S. as a consequence of climate change and increasing human demand for water (Baron et al., 2002; Golladay et al., 2016). This region also contains the highest diversity of freshwater mussels globally (Williams et al., 1993; Haag, 2010). Thus, future mussel mass mortality events are highly likely and we need to both understand their ecological effects and how to mitigate them. Individual mussel species' tolerances to maximum water temperature and minimum dissolved oxygen concentrations vary and are an area of active research (Archambault et al., 2014; Jeffrey et al., 2018). Understanding how mussels acclimate and potentially adapt to increased water temperatures and reduced water availability will be critical to protecting this diverse guild (Galbraith et al., 2012; Gough et al., 2012). However, it is unlikely that the drought sensitive species will rebound to their former abundance without changes in water management and restoration of populations through mussel propagation.

Freshwater mussels are not the only organism threatened by mass mortality events and rivers are not the only ecosystem altered through these events. Our research provides an example of how the loss of an abundant, long-lived organism has cascading and long-term impacts on ecosystems. These impacts are analogous to loss of a forest in terrestrial ecosystems; habitat provision and nutrient sequestration is altered as the community shifts and takes decades to rebound (Ellison et al., 2005; Boyd et al., 2013). The loss of this long-lived organism and the subsequent release of this nutrient pulse has large impacts on stream ecosystems.

Ethics Statement

Fish were handled in accordance with the recommendations of American Fisheries Society's Guidelines for the use of Fishes in Research. The protocol was approved by the University of Oklahoma's Institutional Animal Care and Use Committee.

Author Contributions

TD, CV, and CA designed the mesocosm experiment and TD performed it. CA, CV, and SG provided field data on mussel losses and shell dissolution rates. TD performed analyses and wrote the manuscript with input from all authors.

Funding

Funding to conduct the mesocosm experiment was provided by a University of Oklahoma Biological Station fellowship to TD and NSF DEB-1457542 to CV. Oklahoma mussel surveys were funded by the Oklahoma Biological Survey and an EPA STAR fellowship to CA. Georgia mussel surveys were funded by the Jones Center at Ichauway, Georgia Department of Natural Resources, and The Nature Conservancy.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We appreciate field and laboratory help from K. Ashford, M. Carman, N. Ferreira Rodríguez, J. Hartwell, E. Higgins, J. Lopez, R. Prather, B. van Ee, and M. Winebarger. We thank the Vaughn lab, Atkinson lab, D. Knapp, C. Curry, M. Patten, J. Wesner, and three reviewers for their advice on the manuscript. This paper was completed as part of a dissertation at the University of Oklahoma and is a contribution to the program of the Oklahoma Biological Survey.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2019.00274/full#supplementary-material

References

Allen, D. C., Vaughn, C. C., Kelly, J. F., Cooper, J. T., and Engel, M. H. (2012). Bottom-up biodiversity effects increase resource subsidy flux between ecosystems. Ecology 93, 2165–2174. doi: 10.1890/11-1541.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Archambault, J. M., Cope, W. G., and Kwak, T. J. (2014). Survival and behaviour of juvenile unionid mussels exposed to thermal stress and dewatering in the presence of a sediment temperature gradient. Freshw. Biol. 59, 601–613. doi: 10.1111/fwb.12290