. 2021 Jan 25;118(6):e2008223118. doi:10.1073/pnas.2008223118

Benthic invaders control the phosphorus cycle in the world’s largest freshwater ecosystem

Jiying Li^a,^b,¹,Vadym Ianaiev^a,Audrey Huff^a,John Zalusky^a,Ted Ozersky^a,^c,Sergei Katsev^a,^d,¹

^aLarge Lakes Observatory, University of Minnesota Duluth, Duluth, MN 55812;

^bDepartment of Ocean Science, The Hong Kong University of Science and Technology, Hong Kong, China;

^cDepartment of Biology, University of Minnesota Duluth, Duluth, MN 55812;

^dDepartment of Physics and Astronomy, University of Minnesota Duluth, Duluth, MN 55812

To whom correspondence may be addressed. Email:skatsev@d.umn.edu orlixx0590@d.umn.edu.

Edited by David Strayer, Cary Institute of Ecosystem Studies, Ann Arbor, MI, and accepted by Editorial Board member Mary E. Power November 19, 2020 (received for review April 29, 2020)

Author contributions: J.L. designed research; J.L., V.I., A.H., J.Z., T.O., and S.K. performed research; J.L., V.I., A.H., J.Z., T.O., and S.K. analyzed data; and J.L., T.O., and S.K. wrote the paper.

Series information

From the Cover

Issue date 2021 Feb 9.

Published under thePNAS license.

PMC Copyright notice

PMCID: PMC8017983 PMID:33495360

Significance

The ecosystems of the Laurentian Great Lakes have been dramatically reengineered by invasive bottom-dwelling dreissenid mussels. A key question is whether this biological change has altered whole-ecosystem biogeochemistry, in particular the cycling of phosphorus, a key nutrient that limits biological productivity in freshwater ecosystems. We show that phosphorus cycling in the invaded Great Lakes is now regulated by the population dynamics of a single benthic species, the quagga mussel. This qualitatively changes the responses of the affected lakes to phosphorus inputs from watersheds, complicates predictions, and necessitates a new paradigm for managing these large aquatic ecosystems. Similar changes likely play out in many other dreissenid-invaded lakes across Europe and North America.

Keywords: invasive species, dreissenid mussels, phosphorus cycle, the Great Lakes

Abstract

The productivity of aquatic ecosystems depends on the supply of limiting nutrients. The invasion of the Laurentian Great Lakes, the world’s largest freshwater ecosystem, by dreissenid (zebra and quagga) mussels has dramatically altered the ecology of these lakes. A key open question is how dreissenids affect the cycling of phosphorus (P), the nutrient that limits productivity in the Great Lakes. We show that a single species, the quagga mussel, is now the primary regulator of P cycling in the lower four Great Lakes. By virtue of their enormous biomass, quagga mussels sequester large quantities of P in their tissues and dramatically intensify benthic P exchanges. Mass balance analysis reveals a previously unrecognized sensitivity of the Great Lakes ecosystem, where P availability is now regulated by the dynamics of mussel populations while the role of the external inputs of phosphorus is suppressed. Our results show that a single invasive species can have dramatic consequences for geochemical cycles even in the world’s largest aquatic ecosystems. The ongoing spread of dreissenids across a multitude of lakes in North America and Europe is likely to affect carbon and nutrient cycling in these systems for many decades, with important implications for water quality management.

Invasive species are well known to impact many aspects of ecosystems, including biodiversity, food web structure, and ecosystem functioning (1). The Laurentian Great Lakes, the largest freshwater ecosystem on Earth, serve as a dramatic example of large-scale reorganization of geochemical cycles by an invader. Following the establishment of zebra and quagga (dreissenid) mussels in littoral areas of the Great Lakes in the late 1980s, nutrients and productivity were redirected to the nearshore (2), while pelagic primary productivity declined by as much as 70% (3–7). Having outcompeted zebra mussels, quagga mussels now are abundant in most of the bottom areas in all of the Great Lakes except Lake Superior, often at densities exceeding 10,000 individuals per square meter (8–11). The expansion of quagga mussels coincided with unexplained changes in the abundances and distributions of other benthos (12) and modifications to the structure and phenology of the phytoplankton community (13) and food web structure (14,15). Less attention was given to observations that pelagic concentrations of phosphorus (P), the productivity-limiting nutrient in the Great Lakes, decreased even while external P inputs remained steady (16,17). The dynamics of P have practical importance because regulation of P inputs from the watershed has been the primary tool for managing water quality in the Great Lakes. In particular, reductions in P loadings are credited for the recovery from eutrophic conditions of the 1970s, and a further 40% reduction has been proposed recently to curtail recurring algal blooms in Lake Erie (18).

However, will the dreissenid-colonized lakes respond to further reductions in P input in the same way they did in the past? The unprecedented changes in the Great Lakes induced by dreissenids warrant a reevaluation of P cycling and budgets, which need to account for mussel biomass and modified benthic–pelagic exchanges. The problem extends well beyond the Great Lakes: dreissenids have now been documented in thousands of inland lakes and rivers in North America (19) and Europe (20,21) and may affect freshwater nutrient dynamics on continental scales.

The dynamics of P concentrations in lake water are regulated by the balance of P sources and sinks. These include inputs from the watershed, removal with outflow, and net burial in sediments (16,22). The role of the benthic system is significant. Sediments can recycle a large fraction of the deposited P and resupply it to the water column. Before the dreissenid invasion, this recycled fraction of P sedimentation varied from as little as 10% in Lake Michigan to as much as 60% in Lake Erie (16), contributing 15 to 48% of all (internal and external) phosphorus inputs to the water column (16). The efficiency of sedimentary P recycling also determines the inertia with which P concentrations in the water column respond to external inputs (16). When only a small fraction of the deposited P is recycled, total phosphorus (TP) concentrations respond with time lags of only a few years, even in large lakes (16). Recent TP dynamics in the Great Lakes, however, seems to have been decoupled from external P loadings (17), with change beginning shortly after the dreissenid invasion.

Aquatic consumers can strongly impact the cycling of nutrients (23–26). High abundances of dreissenid mussels in particular can modify the sediment–water exchanges of phosphorus. As epibenthic filter feeders, dreissenids continually remove particulate P from bottom water (27). Most (∼90%) of the ingested phosphorus is remobilized (28): apart from a small portion incorporated into soft tissue and shell, P is either excreted into the water column in dissolved form or egested onto the sediment surface as feces and pseudofeces (27,29,30). The egested P is quickly remineralized to dissolved P via microbial decomposition (31). Mussel biomass P becomes mobilized over longer time scales through decomposition of dead tissues, release of gametes, and dissolution of shells (32). All of these processes modify the natural exchanges of P between sediments and water column, potentially affecting whole-system P mass balance. However, the effects of dreissenids on the sedimentary P fluxes have been evaluated in the Great Lakes only locally (27,29,30), and possibilities of regime-changing tipping points (33) or hysteresis in the geochemical dynamics have not been explored.

Here, we show that the increase of invasive mussel biomass and consequent enhancement of benthic P fluxes have pushed the Great Lakes into a new dynamic regime where P concentrations and fluxes are controlled by mussel physiology and population dynamics, while responses to external phosphorus inputs have become muted. This regime is sensitive to biological perturbations, responds more slowly to external regulation, and presents unique challenges to managing these large ecosystems.

The Mass Balance Approach

We analyzed phosphorus mass balance with a dynamic model (16), which was modified to explicitly consider the cycling of P through mussel populations (Table 1,Materials and Methods, andSI Appendix). In addition to simulating the concentrations of total P in the water column (W), the model tracked P amounts in living (M_L) and dead (M_D) mussel tissues, mussel shells (live and dead shells,M_Sh), and egesta (feces and pseudofeces,M_Eg). Besides the traditional flux components of lake models—P inputs from watershed (I), removal with outflow (O), settling with particles into sediments (F_sed.in), and recycling from sediments (F_sed.out)—the model considered physiological fluxes: P uptake into mussel biomass through filter feeding (F_m.in), P released through mussel excretion (F_m.ex), egestion (E_m), and remineralization of dead biomass and shells (F_D andF_Sh). Physiological flux rates may be affected by the life stages and size distributions of dreissenid mussels. For example, mineralization of gametes and nonviable veliger larvae may release P faster into the water column than mineralization of the equivalent biomass of mature mussel remains, and small mussels have higher mass-specific excretion rates than large mussels (34,35). As age and size distributions of mussels vary across lake environments (10,34), physiological parameters such as water clearance rate (the volume of water filtered through mussel bodies per unit time) and mussel mortality rates were selected individually for each lake (Table 1 andSI Appendix, Table S2). By necessity, these were treated as average quantities over the entire lake, without explicit ties to mussel sizes, ages, or temporal and spatial heterogeneities in distributions. Such nuances, although potentially important, were not simulated for lack of calibration data.

Table 1.

Model parameters

Parameter	Lake Michigan	Lake Huron	Lake Erie	Lake Ontario
Volume, V (km³)	4,900	3,540	480	1,710
Water residence time, τ (y)	62 (76)	21 (76)	2.7 (76)	7.5 (76)
Water flushing rate, γ (y⁻¹) = 1/τ (y⁻¹)	0.0162	0.0476	0.370	0.133
Average depth, Z (m)	85	59	19	86
Water column particle export efficiency, ω^*	0.24	0.30	SI Appendix, Table S2	0.24
Apparent settling velocity, ν_a (m ⋅ y⁻¹)	19 (77)	17 (77)	37 (16)	19 to 29 (16)
Apparent settling rate, η = ν_a/Z (y⁻¹)	0.22	0.29	3.5	0.29
Sediment recycling efficiency,e_s =F_sed.out/F_sed.in^†	0.15	0.15	0.60	0.15
Mussel egestion efficiency, ε_Eg =E_m/F_m.in	0.36 (27)	0.36 (27)	0.36 (27)	0.36 (27)
Mussel excretion efficiency, ε_Ex =F_m.ex/F_m.in	0.54 (27)	0.54 (27)	0.54 (27)	0.54 (27)
Mussel water clearance rate,R_Cl, (m³ ⋅ y^-1 ⋅ g⁻¹)^‡	71 (<2010)	71	SI Appendix, Table S2	33
Mussel water clearance rate,R_Cl, (m³ ⋅ y^-1 ⋅ g⁻¹)^‡	33 (>2010)	71	SI Appendix, Table S2	33
Mussels mortality rate, λ (y⁻¹)^§	0.20	0.20	0.25	0.27
Phosphorus partition ratio into shells, α_Sh^¶	0.35	0.35	0.35	0.35
Reactivity of mussel remains,k_D, (y⁻¹)^#	18	18	18	18
Reactivity of mussel egesta,k_Eg (y⁻¹)	58 (31)	58 (31)	58 (31)	58 (31)
Reactivity of mussel shells,k_Sh (y⁻¹)	0.4 (32,78)	0.4 (32,78)	0.4 (32,78)	0.4 (32,78)

Open in a new tab

Water column particles reach bottom waters and become accessible to mussels via passive settling or water column mixing. Particle export efficiency in Lake Michigan was calculated as the fraction of particles from primary production that become deposited into sediments using carbon sedimentation rate and average primary productivity (36,37). Estimates for Lakes Huron and Ontario were obtained from an empirical relationship with water depth (38). The chosen values of particle export efficiency produce premussel modeled TP that that are consistent with measurements.

^{^†}

Sediment recycling efficiencies for Lakes Michigan, Huron, and Ontario were assumed to be similar to those in well-oxygenated sediments in Lake Superior (22), which is in the typical range for freshwater lakes (39). Higher recycling was used for Lake Erie where sediments are anoxic or with very shallow oxygen penetration (40).

^{^‡}

Mussel water clearance rateR_Cl (m³ ⋅ y^-1 ⋅ g⁻¹) is a fitting parameter. Literature estimates are in the range 18 to 280 m³ ⋅ y⁻¹ ⋅ g SFDW⁻¹ (18 m³ ⋅ y⁻¹ ⋅ g⁻¹ for zebra mussels at low seston concentrations similar to the Great Lakes (41); 44 m³ ⋅ y⁻¹ ⋅ g⁻¹ for quagga mussels in Lake Michigan at 2 to 7 °C (13); and >75 to 280 for quagga mussels in Lake Michigan estimated from egested biogenic silica (29)). Our modeledR_Cl is at the lower end of the range. The populations of larger sized mussels (>10 mm) in Lake Michigan were considerably larger in 2015 compared to 2010, which likely contributed to the continued increase of biomass despite declines in population density (10). We adjusted the water clearance rate accounting for this effect.

^{^§}

The fraction M(t)/M(0)= e^−λt is the probability that a mussel will reach aget. The mean residence time (life span) of a mussel is therefore the integral over ages 0 to infinity, or life span = $\int_{0}^{\infty} e^{- λt} dt = 1 / λ$ . The mortality rate of 0.20 to 0.30 y⁻¹ corresponds to a life span of ∼3 to 5 y. The life span of zebra mussels was estimated to be 4 to 7 y (λ = 0.14 to 0.25 y⁻¹), specific to environments (42). We use an average life span of 4 to 5 y for Lakes Michigan, Huron, and Ontario (λ = 0.2 to 0.30 y⁻¹). The mortality constant should indirectly account for the release of gametes and nonviable veligers, as we do not model this separately. Mortality rates between 2005 and 2010 in Lake Ontario were increased to 1 y⁻¹, taking into account the increases in predation pressure caused by invasive round gobies (11,43).

^{^¶}

The mass ratio of shell:tissue (dried) is about 16.8 (dry tissue accounts for 5.6% of TDW (44); also measured to be ∼6.0% in zebra mussels). Therefore, the shell = 0.94 TDW and SFDW = 0.056 TDW. P composition in soft tissue and shells was measured to be 6.2 mg/g SFDW and 0.15 mg/g shell on average, respectively (SI Appendix, Table S4). Therefore, P composition in shells = 0.15 mg/g shell = 0.15 mg/g × 0.94 TDW = 0.141 mg/g TDW; P composition in tissue = 6.2 mg/g SFDW = 6.2 mg/g × 0.056 TDW = 0.35 mg/g TDW. Therefore, the P partition ratio into the shell is α_Sh = 0.35.

^{^#}

Dead mussel tissues decompose within several weeks. Using 2 wk as half-time,k_D = −ln(1/2)/(14 d) = 0.0495 d⁻¹ = 18 y⁻¹. The model is not sensitive tok_D.

We present the results from both a linear model (Materials and Methods,Table 1, andSI Appendix, section S2), in which mussel growth is determined by food availability but not by the population density or habitat availability, and a nonlinear (logistic) model, which accounts for increased mortality at high population densities (Materials and Methods andSI Appendix, section S3 and Table S2). The linear model is amenable to rigorous mathematical analysis and demonstrates several important insights, which help define some key limits for the water column phosphorus concentration. The dynamics simulated with the linear model are applicable for the initial expansion of mussels, whereas the logistic model better describes the phosphorus dynamics if mussel populations eventually stabilize.

P inventories and fluxes were modeled on multiannual time scales that are sufficiently long so that seasonal and spatial patterns within the lake do not need to be resolved, and P pools can be represented by their annual averages. As the Great Lakes mix vertically every year, these conditions are generally fulfilled over the 2- to 3-y time spans on which the waters are exchanged laterally between basins (17). As our focus was on factors common to all dreissenid-invaded Great Lakes, influences important in individual lakes were not included in the model but are discussed separately. For example, in Lake Erie, hypoxic events induce mussel mortality in the deeper central basin but not in the shallow western basin (45). This effect is taken into account by adjusting the mortality rates in Lake Erie (Materials and Methods andSI Appendix, section S4 and Table S2). In Lake Ontario, mussel mortality induced by upwelling and predation by invasive round gobies (11) were considered explicitly (Table 1). Model results were calibrated to literature data and to the results of our own measurements of benthic fluxes and mussel excretion rates, which were quantified experimentally in Lakes Michigan and Huron over a range of water depths and mussel population densities (Materials and Methods andSI Appendix, Table S3).

Nondreissenid biota were not included in our mass balance analysis. While the abundance of several benthic taxa, including native unionid mussels,Diporeia spp. amphipods, and oligochaete worms, changed following dreissenid establishment (46,47), their preinvasion biomass was substantially lower than that of dreissenid mussels, which now strongly dominate (>90%) total benthic animal biomass in the Great Lakes (7,8,12,47–49). Therefore, changes in P dynamics caused by changes in abundances of nondreissenid benthic consumers are likely small relative to the role of dreissenids.

Results

Comparing P pools and fluxes before and after the invasion (Figs. 1 and2) reveals that dreissenid mussels have become a major agent of P cycling in all four of the invaded Great Lakes. The tissues and shells of quagga mussels now contain nearly as much phosphorus as the entire water columns of the impacted Great Lakes (Figs. 1 and2 andSI Appendix, Table S5). The benthic fluxes of phosphorus scale with mussel biomass (Figs. 1D and2D) and the approximately linear relationship used by the model (Eqs.2–5) falls within the range suggested by the site-specific sediment incubations (Fig. 3 andSI Appendix, Table S6) and the biomass-specific metabolic rates reported in literature (27,29). Critically, mussel-colonized sediments exchange P with the water column orders of magnitude faster than mussel-free sediments (Figs. 1D and2D). In Lake Michigan, for example, dreissenid filter feeding outpaced passive sedimentation as a sink for P around the year 2000 and now exceeds it >10-fold (Fig. 1D). Excretion and egestion now resupply P from sediments into the water column an order of magnitude faster than the preinvasion transport by molecular diffusion and bioirrigation (Fig. 1D). At ∼8 times the rate at which P enters the lake from watershed (Fig. 4), these fluxes now dominate not only the sediment–water exchanges of P but also the P balance for the entire lake (Fig. 4).

Fig. 1. — Temporal dynamics of phosphorus pools and fluxes in Lakes Michigan, Huron, and Ontario. (A) External P inputs (I) (60); (B) TP concentration in the water column during spring overturn. Symbols are measured data from the literature [☐ (6,17); ▽ (7)]. Lines are model results: the dynamic simulation (solid black), the theoretical equilibrium value for a mussel-free system (W_ST; blue dotted), and the asymptotic value after mussel invasion (W_M; red dotted). (C) Phosphorus amounts within the mussel pools, converted to concentration units using the respective lake volumes to simplify comparisons. P inventories in live mussels (red stars) were calculated from measured mussel biomasses (8–11) using the measured 6.2 mg P per gram of SFDW (*SI Appendix*, Table S4). (D) The temporal evolution of benthic fluxes. Shading marks the regions of extrapolation. Note the different scales for they-axes.

Fig. 2. — Temporal dynamics of phosphorus pools and fluxes in the three basins of Lake Erie. (A) External P inputs (I) (18,60,72,75). (B) Total phosphorus (TP) concentration in the water column during spring overturn. Symbols are measured data from the literature (17). Lines are model results: the dynamic simulation (solid black), the theoretical equilibrium value for a mussel-free system (W_ST; blue dotted); and the asymptotic value after mussel invasion (W_M; red dotted). (C) Phosphorus amounts within the mussel pools, converted to concentration units using the respective lake volumes, to simplify comparisons. P inventories in live mussels (red stars) were calculated from measured mussel biomasses (10,45,74) using the measured 6.2 mg P per gram of shell-free dry weight (*SI Appendix*, Table S4). (D) The temporal evolution of benthic fluxes. Shading marks the regions of extrapolation. Note the different scales for they-axes.

Fig. 3. — Relationship between benthic P effluxes and phosphorus amounts contained in mussel tissues. Solid line is the relationship extracted from dynamic mass balance simulations; symbols are site-specific measurements from incubations of mussels (Δ; excretion + egestion) and intact sediment cores (with [open circles]) and without mussels [filled circles]) (seeMaterials and Methods).

Fig. 4. — P budget in Lake Michigan (estimated for year 2020). P pools are shown in parentheses and fluxes are in gigagrams per year. For other lakes, seeSI Appendix, Tables S5 and S7.

The net effect of the uptake of P and its recycling by dreissenids (Figs. 1D and2D) depends on whether the mussel populations are growing or declining. A growing population can rapidly deplete phosphorus from the water column. In Lake Michigan, sequestration into mussel tissues and shells accounted for 20 to 40% of the total benthic P sink since around 2010 (Fig. 1D andSI Appendix, Table S7). Biomass increase may have subsequently decelerated (Figs. 1C and2C) as geographic expansion slowed down and population composition shifted to larger and older mussels (10). Larger mussels have lower mass-specific filtration and growth rates (50). When mussel populations decline, phosphorus from decaying biomass may be quickly resupplied into the lake. In the mid 2000s, mussel population growth reversed temporarily in lakes Ontario and Erie. In Lake Ontario, biomass declined from 2003 to 2008 in shallow (10 to 30 m) areas, coincident with increases in the populations of invasive round gobies, which can exert strong predation pressure on dreissenids (11,43) (Fig. 1). In Lake Erie, mussel biomass fluctuated dramatically, especially in the western and central basins (Fig. 2), likely in response to hypoxia (51,52). Simulations of these effects (Figs. 1 and2) show that mortality can quickly increase the TP concentrations in the water column, while resumed population growth has the opposite effect.

Discussion

Animals can influence nutrient dynamics in aquatic ecosystems by sequestering, remineralizing, or physically moving nutrients between habitats (23–25,53). Whereas the biogeochemical consequences of animal invasions, including those by dreissenids, have received some attention (26,44,54,55), important gaps in understanding the whole-ecosystem biogeochemical effects of invaders remain (56–58). Here, we show how an invasive bivalve has reengineered nutrient dynamics in the world’s largest freshwater ecosystem, with large implications for ecosystem functioning and water quality management.

A New Dynamic Regime.

Quagga mussels have qualitatively changed the phosphorus dynamics in the Great Lakes, in terms of both phosphorus availability in the water column and the time scales on which it fluctuates and adjusts to changes in P inputs. The difference can be illustrated by analyzing the specific factors that determine P availability in the water column. TP concentrations in dreissenid-free lakes reflect mainly the balance between external inputs and P removal to outflow and net sedimentation (SI Appendix, section S1 andEq.15). If conditions remain approximately constant, within several years, TP concentrations approach their theoretical steady state value (W_ST), which can be calculated fromEq.15.

In contrast, the postinvasion P concentrations in the Great Lakes respond strongly to changes in mussel biomass and are continually affected by mussel population dynamics. The dynamics of the mussel populations and phosphorus concentrations are interlinked. One potential feedback involves limitation of dreissenid growth by the availability of phytoplankton food (53). As the abundance of phytoplankton is tied to the availability of the productivity-limiting nutrient (P), mussel growth rates may be assumed approximately proportional to the water column concentration of phosphorus (W) (Eqs.2 and9). Mussel populations then may stabilize upon reaching a low enough concentration of phosphorus (W_M;Eq.17). When the balance between growth and mortality is therefore determined by mussel physiology (Eq.17 in our formulation; see alsoSI Appendix, section S2), the asymptotic concentration after mussel invasionW_M is different from the predreissenid equilibrium TP levelW_ST (Eq.15). WhileW_ST is defined by the physical and biogeochemical characteristics of the water column and sediment, the concentration limitW_M in our model is dependent on mussel physiology and is rather independent of the external P loading (as inEq.17). The system then has the capacity for two qualitatively different dynamics. One reflects the preinvasion regime where fluxes induced by mussels are negligible, and P concentration is stabilized atW_ST by the preinvasion mass balance (Eq.15). Another is controlled by mussel physiology and population dynamics: mussels sequester P until their growth becomes limited, and phosphorus concentration in the water column (W) asymptotically approachesW_M (Eq.17; seeSI Appendix, section S2 for detailed explanation). When mussels proliferate, therefore, biological fluxes begin to dominate P cycling, and the externally controlled regime yields to the regime that is controlled by the dynamics of mussel populations.

Of course, mussel biomass may be constrained by factors other than food supply, such as available sediment space or the effects of high population densities on growth or mortality. Our simulations that consider a logistic increase in mortality (seeMaterials and Methods andSI Appendix, section S3) illustrate that when the population stabilizes (Fig. 5, and as in Lake Erie,Fig. 2), TP concentrations begin approachingW_ST again (Fig. 5) (seeSI Appendix, section S3 for mathematical explanation). In a crucial difference from the mussel-free regime, however, the approach toW_ST takes significantly longer than before the invasion (Fig. 5), as some of the externally supplied P continues to be directed into the biomass. The results indicate that the timing and the trajectory of the approach are strong functions of poorly constrained physiological parameters, such as the degree of nonlinearity in the mussel mortality rate (see alsoSI Appendix, Fig. S4 for a simulated trajectory following a rapid collapse in mussel biomass).

Fig. 5. — Effects of mussel growth dynamics on TP concentrations (for Lake Michigan). (A) The dynamics of TP concentration (W) in the water column. Dashed black: under mussel-free conditions, TP quickly approaches the theoretical steady state valueW_ST (dotted blue) controlled by external loading. Solid black: when mussel populations grow unchecked (linear growth model),W approaches the fixed asymptoteW_M (dotted red). Solid gray: under bounded logistic growth (λ = λ₀[1 + qM_L]; q = 0.02 Gg⁻¹), the asymptoteW_M (dotted purple) itself evolves, eventually approachingW_ST. When mussels are present, TP concentrations approach their respective equilibria much more slowly than in the absence of mussels. (here, the external loading of P and water clearance rate were the same as inFig. 1). (B) TP amount (converted to concentration units using the volume of Lake Michigan) in mussel tissues for the two growth models: the linear model (solid red) and the logistic model (solid purple). Red stars correspond to data calculated from measured mussel biomasses (8–11) using the measured value of 6.2 mg P per gram of SFDW (*SI Appendix*, Table S4).

This analysis allows us to obtain the theoretical limits for the TP trajectories in the Great Lakes, as mass balance dictates that they should remain within the range bounded byW_M andW_ST (seeSI Appendix for details). If quagga mussel populations continue to increase, TP concentrations in Lakes Michigan and Huron may be expected to decrease further, stabilizing in the 2 to 3 and 1 to 2µg/L range, respectively, corresponding to ultraoligotrophic conditions (Fig. 1). In Lake Ontario, TP concentrations would be expected to stabilize around 4µg/L, slightly below their present levels of 6µg/L (Fig. 1). The decline in TP may stop if mussel expansion into the profundal areas slows down or biomass stabilizes for some other reason. In that case, TP concentrations may slowly rebound to theW_ST values that reflect external P loadings: as high as 4µg/L in Lake Michigan, 3µg/L in Lake Huron, and 8µg/L in Lake Ontario (Fig. 1). As mussel biomass no longer increases in Lake Erie, the dreissenid effect on TP concentrations in the lake is weakening (Fig. 2), and the relatively short hydrological residence time (∼3 y) should keep this lake responsive to external P inputs.

Reduced Sensitivity to External P Inputs.

The cycling of large amounts of P through quagga mussel populations transfers the control of P dynamics away from external inputs. Mussel growth (or mortality) can take up (or release) large quantities of P, offsetting or overshadowing the effects of external loading. For example, in Lake Michigan, the entire annual input of P from watershed (4 Gg P) may be absorbed by an increase in mussel biomass of less than 50% (Fig. 4). Unimpeded growth of the population, if realized, would continue depleting TP from the water column even if the external P inputs returned to their 1970s levels (SI Appendix, Fig. S1A). Rapid increases in mussel biomass are not unrealistic, as they have been observed in the Great Lakes as well as in other systems (59). It is not yet clear at what levels the dreissenid populations will stabilize in Lakes Michigan and Huron, nor whether they are likely to experience strong fluctuations, such as the recent ones in Lakes Ontario and Erie (10,11).

One may explore a theoretical possibility of curtailing mussel growth by limiting the external inputs of P to levels where the decreased primary productivity would cause their starvation. Below a certain threshold, the water column TP is expected to drop below the minimal level (W_M) required for mussel growth (SI Appendix, section S2 and Fig. S1B). TP would then be again controlled by external input, approachingW_ST (Eq.15 andSI Appendix, Fig. S1B). For Lakes Michigan and Huron, achieving this is impractical, as the external inputs would need to be cut 60 to 75% to ∼1.4 Gg P/y (SI Appendix, Fig. S2). Atmospheric inputs alone contribute ∼0.3 Gg P/y in Lake Michigan (60), and combined atmospheric and inflow contributions contribute 0.8 Gg P/y to Lake Huron (60). In Lake Ontario, the threshold value is higher, at 2.2 Gg P/y, and theoretically may be realized (SI Appendix, Fig. S2). Even if the external inputs were to drop, however, the legacy P within the mussel biomass would prolong the response to over 20 y (SI Appendix, Fig. S1B). In Lake Erie, nevertheless, reducing external inputs is likely to limit mussel growth over a reasonably short time scale of several years (SI Appendix, Fig. S3). This has practical significance, as the 2019 US–Canada agreement calls for a 40% reduction in TP loadings to the western and central basins of Lake Erie, from the 2008 baseline, to combat algal blooms and hypoxia. In the case of Lake Erie, such a reduction appears to have a high chance of success in controlling TP levels.

The weakened response to external loading complicates the management of invaded lakes, by decreasing the effectiveness of one of the primary tools for regulating their trophic status. Monitoring and control of external inputs, however, remain important, both because the external inputs become important again if mussel populations eventually stabilize (Fig. 5) and because nutrient loading may increase mussel growth (SI Appendix, Fig. S1A (53)). Even if the phosphorus concentrations in the water column remain low, accumulation of P in mussel biomass (SI Appendix, Fig. S1A) creates a high risk of large and unpredictable changes in the system if mussel biomass suddenly declines (e.g., as in Lakes Erie and Ontario,Figs. 1 and2), releasing P into the water column (SI Appendix, Fig. S4).

The New Unpredictable Future.

Our results show that introduction of a single species may shift even large aquatic ecosystems into new dynamic regimes from which they may be slow to return. Quagga mussels have changed P cycling in the Great Lakes by transiently storing large quantities of P in their biomass, drastically modifying the rates of benthic–pelagic exchanges through their biological activities and affecting the long-term properties of sediments by modifying their chemical composition. The phosphorus dynamics in the Great Lakes now depend strongly on the dynamics of invasive mussel populations and less on external P inputs. These changes make this large ecosystem less manageable and portends larger uncertainties. As evident in recent data, large fluctuations in TP may occur for reasons that may be difficult to identify and control (Figs. 1 and2). Meta-analysis of invaded systems indicates a wide range of dreissenid population dynamics (59). In some systems, mussel populations undergo large periodic cycles (35); in others, they stabilize, and yet some experience booms, substantial declines, or irregular fluctuations (20,61–63). These dynamics may be driven by system-specific and often unpredictable ecological characteristics, such as predation (64), current and wave scour (65), hypoxia (51), extreme temperatures (66), changes in food resources, or competition withinDreissena populations or with new invaders (67–69). This poor predictability of mussel populations now complicates projections for whole-ecosystem productivity.

Regulation of phosphorus dynamics by benthic invaders may also lead to broader geochemical effects. As the productivity-limiting nutrient, phosphorus influences the cycles of other key elements such as carbon and nitrogen. Phosphorus scarcity leads to oligotrophication, while mussel activities redirect carbon from water column into the benthic system. Redirection of organic material into the sediments, the hotspots of nitrogen removal (70), may well result in long-term changes in the nitrogen cycle. By controlling geochemical cycling, mussels now control multiple ecosystem processes, including changes in water quality, primary and secondary production, phytoplankton community and phenology, and food web structure (3,6,7,13–15). Similar scenarios likely play out in many other large and small invaded lakes across the continent (44).

Our analysis indicates several ways in which the understanding of the Great Lakes ecosystem can be improved under the new, mussel-dominated, dynamic regime. While monitoring of external phosphorus loads remains important, routine tracking of mussel biomass and distributions (59,71) should be continued and can be improved by increasing the spatial resolution and reporting times of the ongoing programs. Mussel monitoring programs should also continue to track and report changes in mussel demographics, particularly size distributions and changes in abundance of veliger larvae and juvenile mussels. These data are needed to constrain mussel fertility and mortality parameters, which our modeling indicates strongly affect P dynamics. These data will also help determine potential limits to mussel biomass growth due to factors such as substrate availability or predation pressure. Further biological research is needed to constrain the physiological parameters that determine the theoretical limit for the whole-ecosystem phosphorus concentration (W_M) and for which only limited data are currently available (27,44). These parameters include the rates of dreissenid filtration, egestion, excretion, and shell production, the elemental composition of mussel tissues and biodeposits, and the variations of these parameters with environmental conditions. Finally, an aspect that was largely outside of the scope of our present analysis, but which may ultimately determine the new equilibrium in benthic–pelagic P exchanges, is the effect of dreissenids on sediment biogeochemistry. Mussel biodeposits can change the reactivity of sediment organic matter, dense mussel colonies can physically limit the penetration of oxygen into sediments, and extirpation of nativeDiporeia amphipods modifies bioturbation patterns. These changes can strongly affect the recycling of phosphorus in the sediments and in the long term may become a key factor affecting the phosphorus dynamics after dreissenid populations stabilize. Characterizing these unknowns will help predict the future of the Great Lakes ecosystems in this new era of invader-controlled geochemical cycles.

Materials and Methods

The Model.

The mass balance model was adapted from ref.16. As the Great Lakes mix annually, on the multiannual time scale, the water body can be considered well mixed. Accordingly, the model does not resolve spatial variability within the water column (except in the three distinct basins of Lake Erie), nor does it simulate seasonal dynamics: model variables should be considered as annual averages within the entire lake (SI Appendix). Water column phosphorus concentrations are approximated by the TP measured during spring overturn. Mussel mortality rates, water clearance rates (the rate at which mussels pass water through their bodies while filter feeding), and other physiological parameters are also treated as average values for each lake. They are not explicitly dependent on mussel sizes, ages, or morphs. Where appropriate, unit conversions between inventories (grams) and concentrations (grams per cubic meter) were made using the respective lake volumes (Table 1 andSI Appendix).

Changes in P inventories (grams) were tracked in five pools: water column (W), living mussels (M_L), mussel tissue remains (dead mussel tissues;M_D), mussel shells (M_Sh; including shells of live and dead mussels), and mussel egesta (feces and pseudofeces;M_Eg). The temporal dynamics are described by a system of coupled ordinary differential equations:

\frac{d W}{d t} = I (t) - O - F_{sed . in} + F_{sed . out} - F_{m.in} + F_{m.ex} + F_{D} + F_{Eg} + F_{Sh},

[1]

\frac{d M_{L}}{d t} = (1 - α_{Sh}) (F_{m . in} - F_{m . ex} - E_{m}) - λ M_{L},

[2]

\frac{d M_{Sh}}{d t} = α_{S h} (F_{m . in} - F_{m . ex} - E_{m}) - F_{Sh},

[3]

\frac{d M_{D}}{d t} = λ M_{L} - F_{D},

[4]

\frac{d M_{Eg}}{d t} = E_{m} - F_{Eg} .

[5]

Here, the fluxes on the right (grams per year) are external input (I) of P from watershed, P removal with outflow (O), sedimentation (passive particle settling into the sedimentsF_sed.in), flux out of sediment into the water column (F_sed.out), uptake through mussel feeding (F_m.in), mussel excretion of P (F_m.ex), mussel egestion of P (E_m), P release from degradation of mussel remains and egesta (F_D andF_Eg), and dissolution of shells (F_Sh). Mussel mortality is described by a mortality rate constant (λ). Partitioning coefficient (α_Sh) describes the proportion of P incorporated into shells.

External phosphorus inputs,I(t), including point sources, tributary, upstream contributions (from Lake Michigan to Huron via the Straits of Mackinac, Lake Huron to Erie via the St. Clair River, and Lake Erie to Ontario via the Niagara River), and atmospheric inputs, were taken from literature (60,72,73). Outflow is calculated as proportional to nutrient inventory (concentration) in the lake:

O = γ W,

[6]

where γ (year⁻¹) is the inverse of the hydrological residence time. The average sedimentation of phosphorus with passive particle settling is treated similarly to ref.16 and considered proportional to the nutrient’s total amount in the water column:

F_{sed . in} = η W = \frac{ν_{a}}{Z} W,

[7]

whereη is the apparent settling rate, ν_a is the apparent settling velocity, andZ is the average depth. The return phosphorus flux from sediments,F_sed.out, is defined as the fraction of sedimentation flux,F_sed.in:

F_{sed . out} = e_{s} F_{sed . in},

[8]

wheree_s is the sediment recycling efficiency. The effects of mussels on sediment recycling efficiency are not known. For linear simulation, we keep this parameter a constant.

Mussel ingestion (by filter feeding) depends on food availability and is assumed proportional to the concentration of phosphorusW, the productivity-limiting nutrient. The amount of P assimilated per unit time per gram P of mussel biomass is proportional to the water clearance rate by mussels,R_Cl (per year per gram), scaled by the fraction of particles that escape decomposition in the water column and reach the bottom, ω:

F_{m . in} = ω R_{Cl} W M_{L} .

[9]

Water clearance rates were adjusted to individual lake conditions (Table 1 andSI Appendix, section S2). Excretion (of dissolved P) and egestion (production of biodeposits) were described as fixed fractions of the ingestion flux:

F_{m . ex} = ε_{Ex} F_{m . in},

[10]

E_{m} = ε_{Eg} F_{m . in} .

[11]

Release of P from remineralization of dead tissues, egesta, and shells were described using respective kinetic parameters (year⁻¹) for reactivities of mussel remains (k_D), egesta (k_Eg), and shells (k_Sh):

F_{D} = k_{D} M_{D},

[12]

F_{Eg} = k_{Eg} M_{Eg},

[13]

F_{Sh} = k_{Sh} M_{Sh} .

[14]

The values of all model parameters are listed inTable 1 andSI Appendix, Table S2. The model was calibrated to literature data on P concentrations and mussel biomass by adjusting parameters to fit the temporal trends over the past 50 y.

The Asymptotic Steady States.

For given external loadingI, in the absence of mussels, the TP concentration in the water column approaches the steady state value (see ref.16 andSI Appendix, section S1 for details):

W_{ST} = \frac{I}{(γ + η - e_{s} η)} .

[15]

When mussels are present, the water column phosphorus concentration (W) cannot stabilize if the biomass changes becauseEq.1 depends on mussel biomass. A stable biomass would require balancing growth and mortality, which in the model can be expressed (combiningEqs.2 and9–11) as follows:

\frac{d M_{L}}{d t} = (1 - α_{Sh}) ω R_{Cl} (1 - ε_{Ex} - ε_{Eg}) W M_{L} - λ M_{L} = 0 .

[16]

In a linear approximation where the rate constants (ε_Eg, ε_Ex,R_Cl, and λ) are independent of biomassM_L, the balance (achieving zero right-hand side inEq.16) does not define the value ofM_L, which must be determined from other considerations. This balance, however, defines the asymptotic concentration of phosphorus at which the growth stopsSI Appendix, section S2):

W_{M} = \frac{λ}{(1 - α_{Sh}) (1 - ε_{Ex} - ε_{Eg}) ω R_{Cl}} .

[17]

AsW_M is different fromW_ST (Eq.15), achieving the steady state ofEqs.1 and2 simultaneously is impossible. The phosphorus concentrations (W) approachW_ST when the dreissenid-driven fluxes are small, andW_M otherwise. A global steady state may be achieved when the biomass is affected by higher-order factors, such as increases in mortality at higher biomass. We use a logistic model to simulate these effects (SI Appendix, sections S3 and S4).

Analytical methods.

Mussel P content.

Phosphorus content was measured in mussel tissues and shells collected from multiple locations in Lake Michigan (SI Appendix, Table S4) on board the research vessel (R/V)Blue Heron in the summers of 2018 and 2019. Soft tissues were removed from shells and shells were thoroughly cleaned. Tissues and shell samples were dried at 60 °C and ground into fine powder before analyses. P content in the samples was determined spectrophotometrically using a SEAL Analytical AQ400 discrete analyzer following combustion (4 h at 450 °C) and persulfate digestion.

Benthic fluxes.

Benthic fluxes of P were measured in sediment cores collected from several mussel-colonized locations in Lake Michigan (SI Appendix, Table S3) using an Ocean Instruments Multi-Corer. Intact sediment cores of internal diameter 94 mm were incubated at the in situ temperature of 4 °C at 75 to 100% of the in situ oxygen concentrations for up to 6 d, with overlying water being continuously stirred with magnetic stirrers rotating at 33 rpm. Samples of overlying water were withdrawn at regular intervals, filtered through 0.2-µm pore size filter, and analyzed for soluble reactive phosphorus (SRP) by standard molybdenum blue colorimetric method using an ultraviolet-visible spectrometer. The SRP fluxes from sediments obtained from the observed increases in SRP concentrations represent total net effect between P loss to sediment (sorption and uptake) and P release (mussel excretion, degradation of mussel biodeposits [egesta], and sediment P release). To compare experimental data to model results, mussel biomass in the incubation cores were determined by quantifying mussel total dry weight (TDW). P content was calculated using the P content in tissue of 0.35 mg/g TDW (Table 1).

Mussel excretion rates were independently quantified (SI Appendix, Tables S3 and S6) using quagga mussels collected with a Ponar grab. Mussels were cleaned of biofilm and incubated for 3 h in closed containers at the in situ temperature (4 °C) in the dark. SRP analyses were performed by molybdenum blue colorimetric method using a SEAL Analytical AQ400 discrete analyzer. Mussels used in incubations were removed from their shells and dried in a drying oven to determine shell-free dry weight (SFDW). P content was calculated using the P composition in tissue of 6.2 mg P/g SFDW (Table 1 andSI Appendix, Table S4). Mussel egestion was assumed to be 67% of excretion (27).

Supplementary Material

Supplementary File

pnas.2008223118.sapp.pdf^{(847.7KB, pdf)}

Acknowledgments

This work was supported by an NSF grant to T.O. and S.K. (OCE-1737368). We thank the captain and crew of the R/VBlue Heron for assistance in sampling. Comments by Dr. D.L. Strayer and three anonymous reviewers helped to substantially improve this manuscript.

Footnotes

The authors declare no competing interest.

This article is a PNAS Direct Submission. D.S. is a guest editor invited by the Editorial Board.

Seeonline for related content such as Commentaries.

This article contains supporting information online athttps://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2008223118/-/DCSupplemental.

Data Availability.

All study data are included in the article and/orSI Appendix.

References

1.Gallardo B., Clavero M., Sánchez M. I., Vilà M., Global ecological impacts of invasive species in aquatic ecosystems. Glob. Change Biol.22, 151–163 (2016). [DOI] [PubMed] [Google Scholar]
2.Hecky R. E., et al., The nearshore phosphorus shunt: A consequence of ecosystem engineering by dreissenids in the Laurentian Great lakes. Can. J. Fish. Aquat. Sci.61, 1285–1293 (2004). [Google Scholar]
3.Barbiero R. P., Lesht B. M., Warren G. J., Convergence of trophic state and the lower food web in Lakes Huron, Michigan and Superior. J. Great Lakes Res.38, 368–380 (2012). [Google Scholar]
4.Yousef F., Shuchman R., Sayers M., Fahnenstiel G., Henareh A., Water clarity of the Upper Great lakes: Tracking changes between 1998–2012. J. Great Lakes Res.43, 239–247 (2017). [Google Scholar]
5.Evans M. A., Fahnenstiel G., Scavia D., Incidental oligotrophication of North American Great lakes. Environ. Sci. Technol.45, 3297–3303 (2011). [DOI] [PubMed] [Google Scholar]
6.Dove A., Chapra S. C., Long-term trends of nutrients and trophic response variables for the Great Lakes. Limnol. Oceanogr.60, 696–721 (2015). [Google Scholar]
7.Bunnell D. B., et al., Changing ecosystem dynamics in the Laurentian Great lakes: Bottom-up and top-down regulation. Bioscience64, 26–39 (2014). [Google Scholar]
8.Nalepa T. F., Fanslow D. L., Pothoven S. A., A. J.Foley, III, Lang G. A., Long-term trends in benthic macroinvertebrate populations in Lake Huron over the past four decades. J. Great Lakes Res.33, 421–436 (2007). [Google Scholar]
9.Nalepa T. F., Burlakova L. E., Elgin A. K., Karatayev A. Y., "Abundance and biomass of benthic macroinvertebrates in Lake Michigan in 2015, with a summary of temporal trends" (NOAA Technical Memorandum GLERL-175, NOAA Great Lakes Environmental Research Laboratory, Ann Arbor, MI, 2020).
10.Karatayev A. Y., Burlakova L. E., “Lake Erie and Lake Michigan benthos: Cooperative science and monitoring initiative final report” (USGS-GLRI G14AC00263. Great Lakes Center, Buffalo State, Buffalo, NY 2017).
11.Birkett K., Lozano S. J., Rudstam L. G., Long-term trends in Lake Ontario’s benthic macroinvertebrate community from 1994–2008. Aquat. Ecosyst. Health Manage.18, 76–88 (2015). [Google Scholar]
12.Burlakova L. E., et al., The benthic community of the Laurentian Great lakes: Analysis of spatial gradients and temporal trends from 1998-2014. J. Great Lakes Res.44, 600–617 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
13.Vanderploeg H. A., Liebig J. R., Nalepa T. F., Fahnenstiel G. L., Pothoven S. A., Dreissena and the disappearance of the spring phytoplankton bloom in Lake Michigan. J. Great Lakes Res.36, 50–59 (2010). [Google Scholar]
14.Fera S. A., Rennie M. D., Dunlop E. S., Broad shifts in the resource use of a commercially harvested fish following the invasion of dreissenid mussels. Ecology98, 1681–1692 (2017). [DOI] [PubMed] [Google Scholar]
15.Madenjian C. P., et al., Changes in the lake Michigan food web following dreissenid mussel invasions: A synthesis. J. Great Lakes Res.41, 217–231 (2015). [Google Scholar]
16.Katsev S., When large lakes respond fast: A parsimonious model for phosphorus dynamics. J. Great Lakes Res.43, 199–204 (2017). [Google Scholar]
17.Chapra S. C., Dolan D. M., Great lakes total phosphorus revisited: 2. Mass balance modeling. J. Great Lakes Res.38, 741–754 (2012). [Google Scholar]
18.Great Lakes Water Quality Agreement Nutrient Annex Subcommittee , "Lake Erie binational phosphorus reduction strategy" (Great Lakes Water Quality Agreement Nutrient Annex Subcommittee, 2019). [Google Scholar]
19.Benson A. J., “Chronological history of zebra and quagga mussels (Dreissenidae) in North America, 1988-2010” in Quagga and Zebra Mussels: Biology, Impacts, and Control, Nalepa T. F., Schloesser D. W., Eds. (CRC Press, Taylor & Francis, ed. 2, 2014), pp. 9–32. [Google Scholar]
20.Nalepa T. F., et al., Eds., “Spread of the quagga mussel, Dreissena rostriformis bugensis in Western Europe” in Quagga and Zebra Mussels: Biology, Impacts, and Control (CRC Press, Taylor & Francis, ed. 2, 2014), pp. 83–92. [Google Scholar]
21.Orlova M. I., “Origin and spread of quagga mussels (Dreissena rostriformis bugensis) in Eastern Europe with note on size structure of populations” in Quagga and Zebra Mussels: Biology, Impacts and Control, Nalepa T. F., Schloesser D. W., Eds. (CRC Press, Taylor & Francis, ed. 2, 2014), pp. 93–102. [Google Scholar]
22.Li J., Zhang Y., Katsev S., Phosphorus recycling in deeply oxygenated sediments in Lake Superior controlled by organic matter mineralization. Limnol. Oceanogr.63, 1372–1385 (2018). [Google Scholar]
23.Kitchell J. F., et al., Consumer regulation of nutrient cycling. Bioscience29, 28–34 (1979). [Google Scholar]
24.Atkinson C. L., Capps K. A., Rugenski A. T., Vanni M. J., Consumer-driven nutrient dynamics in freshwater ecosystems: From individuals to ecosystems. Biol. Rev. Camb. Philos. Soc.92, 2003–2023 (2017). [DOI] [PubMed] [Google Scholar]
25.Vanni M. J., Nutrient cycling by animals in freshwater ecosystems. Annu. Rev. Ecol. Syst.33, 341–370 (2002). [Google Scholar]
26.Capps K. A., Atkinson C. L., Rugenski A. T., Implications of species addition and decline for nutrient dynamics in fresh waters. Freshw. Sci.34, 485–496 (2015). [Google Scholar]
27.Mosley C., Bootsma H., Phosphorus recycling by profunda quagga mussels (Dreissena rostriformis bugensis) in Lake Michigan. J. Great Lakes Res.41, 38–48 (2015). [Google Scholar]
28.Stoeckmann A. M., Garton D. W., A seasonal energy budget for zebra mussels (Dreissena polymorpha) in western Lake Erie. Can. J. Fish. Aquat. Sci.54, 2743–2751 (1997). [Google Scholar]
29.Mosley C., “Phosphorus recycling by profunda quagga mussels in Lake Michigan,” MS thesis, University of Wisconsin Milwaukee, Milwaukee, WI (2014).
30.Arnott D. L., Vanni M. J., Nitrogen and phosphorus recycling by the zebra mussel (Dreissena polymorpha) in the western basin of Lake Erie. Can. J. Fish. Aquat. Sci.53, 646–659 (1996). [Google Scholar]
31.Giles H., Pilditch C. A., Effects of mussel (Perna canaliculus) biodeposit decomposition on benthic respiration and nutrient fluxes. Mar. Biol.150, 261–271 (2006). [Google Scholar]
32.Strayer D. L., Malcom H. M., Shell decay rates of native and alien freshwater bivalves and implications for habitat engineering. Freshw. Biol.52, 1611–1617 (2007). [Google Scholar]
33.Strayer D. L., Pace M. L., Caraco N. F., Cole J. J., Findlay S. E. G., Hydrology and grazing jointly control a large-river food web. Ecology89, 12–18 (2008). [DOI] [PubMed] [Google Scholar]
34.Baldwin B. S., et al., Erratum: Comparative growth and feeding in zebra and quagga mussels (Dreissena polymorpha and Dreissena bugensis): Implications for North American lakes. Can. J. Fish. Aquat. Sci.60, 1432 (2003). [Google Scholar]
35.Strayer D. L., et al., Long‐term variability and density dependence in Hudson River Dreissena populations. Freshw. Biol.65, 474–489 (2020). [Google Scholar]
36.Thomsen U., Thamdrup B., Stahl D. A., Canfield D. E., Pathways of organic carbon oxidation in a deep lacustrine sediment, Lake Michigan. Limnol. Oceanogr.49, 2046–2057 (2004). [Google Scholar]
37.Shuchman R. A., Sayers M., Fahnenstiel G. L., A model for determining satellite-derived primary productivity estimates for Lake Michigan. J. Great Lakes Res.39, 46–54 (2013). [Google Scholar]
38.Li J., “Sediment diagenesis in large lakes Superior and Malawi, geochemical cycles and budgets and comparisons to marine sediments,” PhD thesis, University of Minnesota, MI (2014).
39.Caraco N., Cole J., Likens G. E., A comparison of phosphorus immobilization in sediments of freshwater and coastal marine systems. Biogeochemistry9, 277–290 (1990). [Google Scholar]
40.Matisoff G., et al., Internal loading of phosphorus in western Lake Erie. J. Great Lakes Res.42, 775–788 (2016). [Google Scholar]
41.Lei J., Payne B. S., Wang S. Y., Filtration dynamics of the zebra mussel, Dreissena polymorpha. Can. J. Fish. Aquat. Sci.53, 29–37 (2011). [Google Scholar]
42.McMahon R. F., Evolutionary and physiological adaptations of aquatic invasive animals: r selection versus resistance. Can. J. Fish. Aquat. Sci.59, 1235–1244 (2002). [Google Scholar]
43.Naddafi R., Rudstam L. G., Does differential predation explain the replacement of zebra by quagga mussels?Freshw. Sci.33, 895–903 (2014). [Google Scholar]
44.Williamson F., Ozersky T., Lake characteristics, population properties and invasion history determine impact of invasive bivalves on Lake Nutrient dynamics. Ecosystems (N. Y.), 10.1007/s10021-019-00371-z (2019). [DOI] [Google Scholar]
45.Karatayev A. Y., et al., Twenty five years of changes in Dreissena spp. populations in Lake Erie. J. Great Lakes Res.40, 550–559 (2014). [Google Scholar]
46.Crail T. D., Krebs R. A., Zanatta D. T., Unionid mussels from nearshore zones of Lake Erie. J. Great Lakes Res.37, 199–202 (2011). [Google Scholar]
47.Nalepa T. F., Fanslow D. L., Lang G. A., Transformation of the offshore benthic community in Lake Michigan: Recent shift from the native amphipod Diporeia spp. to the invasive mussel Dreissena rostriformis bugensis. Freshw. Biol.54, 466–479 (2009). [Google Scholar]
48.Dermott R., Sudden disappearance of the amphipod Diporeia from eastern Lake Ontario, 1993-1995. J. Great Lakes Res.27, 423–433 (2001). [Google Scholar]
49.Nalepa T. F., Manny B. A., Roth J. C., Mozley S. C., Schloesser D. W., Long-term decline in freshwater mussels (bivalvia: Unionidae) of the western basin of Lake Erie. J. Great Lakes Res.17, 214–219 (1991). [Google Scholar]
50.Schneider D. W., A bioenergetics model of zebra mussel, Dreissena polymorpha, growth in the Great Lakes. Can. J. Fish. Aquat. Sci.49, 1406–1416 (2008). [Google Scholar]
51.Karatayev A. Y., et al., Biomonitoring using invasive species in a large lake: Dreissena distribution maps hypoxic zones. J. Great Lakes Res.44, 639–649 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
52.Zhou Y., Michalak A. M., Beletsky D., Rao Y. R., Richards R. P., Record-breaking Lake Erie hypoxia during 2012 drought. Environ. Sci. Technol.49, 800–807 (2015). [DOI] [PubMed] [Google Scholar]
53.Strayer D. L., Understanding how nutrient cycles and freshwater mussels (Unionoida) affect one another. Hydrobiologia735, 277–292 (2014). [Google Scholar]
54.Strayer D. L., Caraco N. F., Cole J. J., Findlay S., Pace M. L., Transformation of freshwater ecosystems by bivalves: A case study of zebra mussels in the Hudson River. Bioscience49, 19–27 (1999). [Google Scholar]
55.Hall R. O., Tank J. L., Dybdahl M. F., Exotic snails dominate nitrogen and carbon cycling in a highly productive stream. Front. Ecol. Environ.1, 407–411 (2003). [Google Scholar]
56.Ehrenfeld J. G., Ecosystem consequences of biological invasions. Annu. Rev. Ecol. Evol. Syst.41, 59–80 (2010). [Google Scholar]
57.Ehrenfeld J. G., Effects of exotic plant invasions on soil nutrient cycling processes. Ecosystems (N. Y.)6, 503–523 (2003). [Google Scholar]
58.Norkko J., et al., A welcome can of worms? Hypoxia mitigation by an invasive species. Glob. Change Biol.18, 422–434 (2012). [Google Scholar]
59.Strayer D. L., et al., Long‐term population dynamics of dreissenid mussels (Dreissena polymorpha and D. rostriformis): A cross‐system analysis. Ecosphere10, e02701 (2019). [Google Scholar]
60.Dolan D. M., Chapra S. C., Great lakes total phosphorus revisited: 1. Loading analysis and update (1994-2008). J. Great Lakes Res.38, 730–740 (2012). [Google Scholar]
61.Burla H., Ribi G., Density variation of the zebra mussel Dreissena polymorpha in Lake Zurich, from 1976 to 1988. Aquat. Sci.60, 145–156 (1998). [Google Scholar]
62.Burlakova L. E., Karatayev A. Y., Padilla D. K., Changes in the distribution and abundance of Dreissena polymorpha within lakes through time. Hydrobiologia571, 133–146 (2006). [Google Scholar]
63.Hetherington A. L., et al., Invader invaded: Population dynamics of zebra mussels (Dreissena polymorpha) and quagga mussels (Dreissena rostriformis bugensis) in polymictic Oneida lake, NY, USA (1992–2013). Biol. Invasions21, 1529–1544 (2019). [Google Scholar]
64.Carlsson N. O. L., Bustamante H., Strayer D. L., Pace M. L., Biotic resistance on the increase: Native predators structure invasive zebra mussel populations. Freshw. Biol.56, 1630–1637 (2011). [Google Scholar]
65.MacIsaac H. J., Population structure of an introduced species (Dreissena polymorpha) along a wave-swept disturbance gradient. Oecologia105, 484–492 (1996). [DOI] [PubMed] [Google Scholar]
66.Boeckman C. J., Bidwell J. R., “Density, growth, and reproduction of zebra mussels (Dreissena polymorpha) in two Oklahoma reservoirs” in Quagga and Zebra Mussels: Biology, Impacts and Control, Nalepa T. F., Schloesser D. W., Eds. (CRC Press, Taylor & Francis, 2015), pp. 369–382. [Google Scholar]
67.Strayer D. L., Malcom H. M., Long-term demography of a zebra mussel (Dreissena polymorpha) population. Freshw. Biol.51, 117–130 (2006). [Google Scholar]
68.Vanderploeg H. A., Johengen T. H., Liebig J. R., Feedback between zebra mussel selective feeding and algal composition affects mussel condition: Did the regime changer pay a price for its success?Freshw. Biol.54, 47–63 (2009). [Google Scholar]
69.van der Velde G., Paffen B. G. P., van den Brink F. W. B., bij de Vaate A., Jenner H. A., Decline of zebra mussel populations in the Rhine - competition between two mass invaders (Dreissena polymorpha and Corophium curvispinum). Naturwissenschaften81, 32–34 (1994). [Google Scholar]
70.Li J., Katsev S., Nitrogen cycling in deeply oxygenated sediments: Results in Lake Superior and implications for marine sediments. Limnol. Oceanogr.59, 465–481 (2014). [Google Scholar]
71.Pergl J., et al., Need for routine tracking of biological invasions. Conserv. Biol.34, 1311–1314 (2020). [DOI] [PubMed] [Google Scholar]
72.Lesht B. M., Fontaine T. D., Dolan D. M., Great Lakes total phosphorus model: Post audit and regionalized sensitivity analysis. J. Great Lakes Res.17, 3–17 (1991). [Google Scholar]
73.Barbiero R. P., Tuchman M. L., Warren G. J., Rockwell D. C., Evidence of recovery from phosphorus enrichment in Lake Michigan. Can. J. Fish. Aquat. Sci.59, 1639–1647 (2002). [Google Scholar]
74.Burlakova L. E., Karatayev A. Y., Pennuto C., Mayer C., Changes in Lake Erie benthos over the last 50 years: Historical perspectives, current status, and main drivers. J. Great Lakes Res.40, 560–573 (2014). [Google Scholar]
75.Maccoux M. J., Dove A., Backus S. M., Dolan D. M., Total and soluble reactive phosphorus loadings to Lake Erie: A detailed accounting by year, basin, country, and tributary. J. Great Lakes Res.42, 1151–1165 (2016). [Google Scholar]
76.Quinn F. H., Hydraulic residence times for the Laurentian Great lakes. J. Great Lakes Res.18, 22–28 (1992). [Google Scholar]
77.Durham R. W., Joshi S. R., Recent sedimentation rates, 210Pb fluxes, and particles settling velocities in Lake Huron, Laurentian Great Lakes. Chem. Geol.31, 53–66 (1980). [Google Scholar]
78.Ozersky T., Evans D. O., Ginn B. K., Invasive mussels modify the cycling, storage and distribution of nutrients and carbon in a large lake. Freshw. Biol.60, 827–843 (2015). [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Supplementary File

pnas.2008223118.sapp.pdf^{(847.7KB, pdf)}

Data Availability Statement

All study data are included in the article and/orSI Appendix.

Movatterモバイル変換

PERMALINK

Benthic invaders control the phosphorus cycle in the world’s largest freshwater ecosystem

Jiying Li

Vadym Ianaiev

Audrey Huff

John Zalusky

Ted Ozersky

Sergei Katsev