Food depletion regulates the demography of invasive dreissenid mussels in a stratified lake

Lake stratification produces sharp gradients in temperature and pelagic resources which have cascading effects on the traits of aquatic populations, including invasive species and their ecosystem impacts. We study the consequences of such common environmental gradients on the demography of quagga mussels, one of the world's most aggressive invasive species. Coupling a series of in situ experiments with a biophysical model of the pelagic community, we quantify mussel growth and recruitment in littoral vs. profundal benthic habitats of eastern Lake Erie. We found that both severe food depletion and cold temperatures in the hypolimnion during summer stratification cause mussels to grow twice as slowly and currently inhibit recruitment in profundal compared to littoral habitats. Together with the high biomass and large mean mussel size found in long‐term monitoring surveys in the profundal habitats, our results imply that mussels successfully colonizing profundal habitats have relatively long lifespans with low growth rates, and therefore lower productivity, compared to mussels in shallow areas. Consequently, the bulk of dreissenid biomass in Lake Erie and other Great Lakes found in the vast but resource‐poor profundal habitats may have very limited impacts on epilimnetic communities compared to mussels in littoral areas. By contrast, mussels in profundal habitats strongly reduce food availability to hypolimnetic communities throughout the year. Thus, our results explain how the ecosystem impacts of sessile freshwater invaders are likely to vary among habitats and lakes depending on the relative size of profundal populations and the extent of stratification.

The physiological activities of individuals both determine and reflect how species interact with and impact their ecosystem. For example, high rates of reproduction, growth, and mortality within a species correspond to a strong per capita energetic role in the ecosystem. Differences in environmental conditions or community composition across habitats can therefore produce variation in both the demographic rates of a species and its ecological role (J onasson 2004;Benton et al. 2006). This variation can be especially pronounced in widespread species with strong capacity for phenotypic plasticity, such as invasive species (Davidson et al. 2011). Consequently, ignoring demographic variation among habitats by focusing detailed population studies on specific habitats can produce a strong, unrecognized bias in the estimated ecological role of a species. The importance of this bias is especially acute for invaders rapidly expanding their distribution range, as the impacts of such species in novel ecosystems can be simultaneously strong and poorly resolved.
Thermal stratification is one particularly common source of habitat heterogeneity in aquatic ecosystems. In lakes with strong stratification during productive summer seasons, species are distributed along environmental gradients from the warm, well-mixed littoral zone, where the surface mixedlayer intersects the bottom, to the deep and cold profundal zone (located below the compensation depth where insufficient light penetrates to support photosynthesis), where bottom habitats and hypolimnetic waters are isolated by the thermocline from the surface waters for much of the growing season. Benthic invertebrates in these ecosystems face particularly large depth-dependent gradients in food resources and temperature that strongly affect their reproduction, growth, and productivity (Brinkhurst 1974;Wetzel 1983;J onasson 2004). Concomitantly, these gradients also govern the distribution of benthic invertebrates and limit most species to a particular depth zone (Brinkhurst 1974;Karatayev et al. 2013).
The byssate bivalve Dreissena rostriformis bugensis (quagga mussel) is a rapidly spreading freshwater invader which creates high densities in a wide range of habitat types, ranging from the very shallow littoral areas to the deepest parts of the lakes (reviewed in Karatayev et al. 1998Karatayev et al. , 2011aKaratayev et al. , 2015. Dreissena species are considered the most aggressive freshwater invaders in the Northern hemisphere (reviewed in Karatayev et al. 2002). They have invaded thousands of waterbodies and typically dominate the biomass of benthic communities (> 90%; Karatayev et al. 1997Karatayev et al. , 2002Karatayev et al. , 2015Higgins and Vander Zanden 2010). Through their feeding and filtering activities, dreissenids alter both ecosystem structure by physically transforming bottom habitats and large-scale ecosystem processes, especially trophic interactions and the availability of food for both pelagic and benthic species (reviewed in Karatayev et al. 2002Karatayev et al. , 2015Keller et al. 2007;Higgins and Vander Zanden 2010). However, due to logistical constraints the demographic parameters (fecundity, growth, longevity) of Dreissena spp. and their ecological effects are typically studied in littoral habitats (i.e., bottom areas of the littoral zone), and predominantly focus on zebra mussels (D. polymorpha). Much less work has been conducted on quagga mussels in profundal habitats (i.e., bottom areas of the profundal zone), where populations appear to grow slowly and show little turnover Karatayev et al. 2018); moreover, even fewer studies explicitly compare populations in littoral and profundal habitats within the same waterbody (reviewed in Karatayev et al. 2014aKaratayev et al. , 2015Nalepa and Schloesser 2014). Consequently, our capacity to understand and predict the ecosystem impacts of quagga mussels is severely limited outside of shallow lakes.
In this article, we quantify the differences in quagga mussel demography among profundal and littoral zones of Lake Erie, and identify the environmental factors driving these differences. This system has the longest history of Dreissena invasion in North America: here, quagga mussels were first detected in 1989 and soon colonized all areas of the lake (Dermott and Munawar 1993;Mills et al. 1993). More recent monitoring has revealed much lower recruitment levels (indicated by the presence of few < 5 mm long mussels) and turnover in population structure in profundal habitats compared to littoral areas (Fig. 1, Karatayev et al. 2018), yet the causes, extent, and ecological implications of such demographic differences are unknown. Here, we couple a series of in situ experiments with a biophysical model of the pelagic community to quantify mussel growth and recruitment rates in littoral vs. profundal habitats of eastern Lake Erie and identify the contribution of environmental factors to observed differences in mussel demography. We also analyze near-bottom turbidity, a surrogate for Dreissena food abundance, from a long-term monitoring program to investigate whether mussel food resources were more abundant early in the invasion (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) than later in the invasion (2007)(2008)(2009)(2010)(2011)(2012)(2013). We compare these food conditions with mussel length-frequency distribution at different stages of invasion to determine whether food depletion was reflected in the mussels' growth and recruitment. Finally, we discuss the implications of these differences in quagga mussel demography for the ecological impacts of this species on the ecosystems of large lakes.

Fig. 1.
Length-frequency distribution of quagga mussels from littoral (< 20 m) and benthic profundal (> 40 m) zones in Lake Erie. Data from Karatayev et al. (2014a), and authors unpublished data, including 2014 Lake Erie CSMI study and data from the four EPA GLNPO Long-term monitoring stations in eastern basin (station ER09, 49 m depth; station ER15M, 63 m depth; station ER63, 46 m depth; station ER93b, 42 m depth) and a littoral station in the central basin (ER93b, 16 m depth) for 2012-2015.

Field experiments
We first compared quagga mussel growth rates in the shallow and deep-water habitats of the eastern basin of Lake Erie in 2015, followed by a second study in 2016 to determine the drivers of observed growth and recruitment gradients across depths in the water column at deep-water locations. This area of Lake Erie develops stable stratification during summer with a large hypolimnion which never goes hypoxic (Karatayev et al. 2018). In both years, mussels were kept in cages made of aluminum 9 3 12 3 9 cm with 2.5 mm holes in each side (Fig. 2). Cages were anchored with two cinderblocks connected by 50 m line (1 cm in diameter). The line was used to locate and remove cages with a hook dragged behind the boat. After retrieving cages at the end of each experiment, the line, cages, floats, and cinderblocks were examined for Dreissena settlement. Young-of-the-year mussels were collected at each site (in 2015) and at each depth (in 2016), the area sampled was recorded, and all collected Dreissena were identified to species, counted, and measured using a digital caliper (0.01 mm). All shell length measurements were rounded to the nearest 0.1 mm. Temperature loggers (HoboV R Onset UA-001-64) were attached to each cage and were set to record temperature every 30 min. To measure individual growth increments, all mussels used in experiments were individually numbered with plastic bee tags (3 mm in diameter) glued to the shell. To detect sizedependent growth rates, we placed 45-49 mussels uniformly spanning a range of initial lengths (8-20 mm) in each cage, with average initial lengths being the same across all cages (Table 1). To quantify the relative magnitude of mussel recruitment across depths in 2016, we sampled newly settled Dreissena uniformly across all depths in the water column and at the bottom by collecting mussels from floats, cages, and polypropylene ropes of our experimental setup, adjusting for substrate area to estimate recruit areal density.
To compare demography in littoral vs. profundal habitats, in 2015 we kept mussels in cages at sites with different depths (littoral habitats: 13.7-14.8 m depth vs. profundal habitats: 48.2-54.3 m depth) for 205 d, from May 15 to December 7 (Fig. 3). We installed two replicate cages at each of the three sites in each habitat. Because substrata were different at the shallow (cobble and sand) and deep (silt) sites, cages were suspended at about 70 cm above the bottom surface with floats to avoid siltation.
Given that both food and temperatures differ strongly between littoral and profundal habitats, in 2016 we examined growth and recruitment along a depth gradient at profundal sites in order to determine the influence of each driver on mussel demography. For this, we placed cages on the bottom and at different depths above the bottom at three profundal sites in eastern basin of Lake Erie for 182 d (09 May to 07 November; Table 1; Fig. 3). At each site, cages were attached to a 36 m line anchored with a cinderblock and a float at the top to keep the chain of cages vertical in the water column. Our uppermost cages were placed 25 m and 35 m above the bottom to measure growth near the thermocline, where gradients in temperature and food availability typically diverge (e.g., due to stratification and the presence of the deep chlorophyll layer), and the role of each factor in determining mussel demography could be resolved. Because the hypolimnion in the eastern basin of Lake Erie starts at about 20 m (Beletsky et al. 2013) and compensation depth at about 25 m (LakeAccess 2006) below the surface, our 35 m above the bottom (16-20 m below surface) cages were located in the epilimnion above the thermocline, while all other cages at the deep water sites were located below the thermocline in the bottom and hypolimnetic water of the profundal zone. Additionally, our 2015 study may have overestimated growth rates in profundal habitats by placing cages above the benthic boundary layer (0.7 m above the bottom). Food resources are likely depleted within the benthic boundary layer due to mussel filtration. Such features are known to limit growth in both freshwater and marine mussels (Fr echette and Bourget 1985;Peterson and Beal 1989;Yu and Culver 1999;Karatayev et al. 2006), including in the Great Lakes (Malkin et al. 2012). To resolve any boundary layer effects on growth and recruitment, we also placed cages directly on the bottom (i.e., on the cinder block anchor) and 0.7 m above the bottom. Given very similar growth in duplicate cages at each site in 2015, we did not duplicate cages at each site and depth, retaining replication only among sites.

Analysis of drivers of demography
Across depths, we evaluated degree-days and the relative phytoplankton concentration throughout our experiment ("relative food" hereafter; estimated using the hydrodynamic model-see Dynamics in mussel food resources) as predictors of demography. In particular, we calculated degree days from  temperature logger data while assuming that the biological zero (i.e., temperature below which growth stops) did not occur in our system; this is motivated by the fact that quagga mussels persist throughout the deep profundal areas of Lake Michigan, where year-round temperatures are consistently below the lowest values measured in our study . For all analyses, we normalized each variable to its mean and standard deviation, and transformed recruit density as log(x 1 1). At the depths of our 2016 experiments, relative food and degree days were strongly correlated (0.98 Pearson correlation). Because mussel recruitment at our sites was absent in recent years but occurred earlier in the invasion (when profundal habitat temperatures were similar but chlorophyll availability was higher), we used forward stepwise Ridge regression in STA-TISTICA (TIBCO Software Inc. 2017) to test whether relative food was a stronger predictor of recruitment than degree days. Given that both temperature and food regulate mussel growth (MacIsaac 1994;Karatayev et al. 2006Karatayev et al. , 2011b, we incorporated results from a prior experiment of temperature effects on mussel growth (Karatayev et al. 2011b) to determine the added effect of the food gradient in our study. This lab experiment used water from eastern Lake Erie in a flowthrough system to estimate a relatively precise effect of temperature on mussel growth rates over the range of temperatures observed in our 2015-2016 study and under constant, realistic food conditions. To account for the fact that mussel growth is strongly size-dependent and that initial mussel sizes differed substantially in our study vs. the 2011b temperature experiment, we first estimated the von Bertalanffy growth rate (K) in every treatment of each study. For both studies, we estimated K as the slope of the linear regression between mussel growth vs. initial size. However, in estimating K for 2011b data, which used only a narrow range of initial sizes, we fixed mussel size at which growth stops (L 1 ) at the average value estimated for our 2015-2016 data. This L 1 value appears to be robust as L 1 estimates were not significantly different among our 2015-2016 experiments and generally consistent with mussel sizes observed in the field (22-28 mm; Fig. 1; Table 1). Finally, we used least squares regression to estimate the mean and variance of the (linear) effect of degree days on growth rates K in the 2011b experiment.
To fit a Bayesian linear model predicting K from degree days, relative food levels, and near-bottom food depletion, we used the rethinking package in R (McElreath 2015; R Core Team 2017). Under strong correlation of two predictors (relative food and degree days), the estimated effect of relative food would not be confounded when the effect of degree days is known exactly (i.e., by simply subtracting degree day effects from observed growth). In our case, the effect of degree days estimated above was relatively precise and incorporated as a strongly informative prior (0.53 6 0.06, mean 6 standard deviation of the Gaussian distribution). We used broad, uninformative priors for the effects of relative food and a boundary layer (a factor denoting whether or not the sample was situated directly on the bottom) bounded by 22 and 2, and likewise for standard deviation of the model residuals (bounded by 0 and 2). Finally, we checked to ensure that the presence of food effects predicted by our model was not sensitive to uncertainty in the effect of degree days and noise in the data (Supporting Information Fig. S1). We include all primary data and code used in throughout the analysis of 2015 and 2016 mussel demography (Supporting Information S3).

Long-term monitoring data
To resolve the spatiotemporal distribution of mussel food resources, we supplemented our analysis with monitoring data to evaluate long-term food resources of Dreissena in Lake Erie as well as a hydrodynamic model of the vertical and seasonal phytoplankton distribution during our study (see next sub-section). To understand the long-term absence of mussel recruitment in profundal habitats, we also analyzed near bottom turbidity specifically to determine whether mussels' food resources early in the invasion (before 2007) were better than later in the invasion (after 2007;Fig. 4). Turbidity measures the concentration of suspended solids in the water, which we assumed to be a surrogate for Dreissena food abundance (mostly detritus, phytoplankton, and bacteria; reviewed in Karatayev et al. 2005). Turbidity data used in this analysis were measured in 1996-2013 at the surface (1 m below surface) and near the bottom (1 m above the bottom) at four eastern basin stations (depths range 32-62 m) of Lake Erie twice a year (during spring and summer surveys) aboard R/V Lake Guardian as part of U.S. Environmental Protection Agency's (EPA) Great Lakes National Program Office (GLNPO) Biology Monitoring Program (Fig. 4). To standardize sampling time across years we used data collected during spring isothermal period (April) and summer stratified period (August) cruises only. At each station samples were collected using a Rosette sampling device (SeaBird Electronics Carousel Water Sampler) deployed from the ship. Measurements were recorded in the on-ship laboratory just subsequent to sampling using a nephelometer, and reported in Nephelometric Turbidity Units (NTU) (Barbiero and Tuchman 2004). Dreissena density and size-frequency distribution was analyzed from EPA GLNPO Biology Monitoring Program stations sampled during summer cruses during 2003-2015.

Biophysical model
To illustrate how vertical mixing, stratification, and mussel grazing influence the dynamics of food availability to mussels at littoral vs. profundal habitats in 2015 and to estimate the relative food availability throughout the water column in 2016, we ran simulations using a 1-dimensional (vertical column) version of the biophysical model described by Rowe et al. (2017). Rowe et al. (2015a) showed that 1dimensional biophysical model simulations can reproduce seasonal patterns of phytoplankton abundance observed in the Great Lakes, including the spring bloom during the late isothermal and early stratification periods, the deep chlorophyll maximum during stratification, and secondary maxima in surface and bottom chlorophyll associated with fall turnover. The biophysical model represented a simplified, phosphorus (P)-limited lower food web, consisting of five functional groups: phytoplankton, zooplankton, dissolved P, detritus, and mussel biomass. Temperature-dependent mussel clearance rate was applied based on Vanderploeg et al. (2010). Additional details are provided by Rowe et al. (2017). Model simulations were initialized on 01 January 2015 or 2016 with total phosphorus set to 6 lgPL 21 , which was the mean value from U.S. EPA survey data in the eastern basin of Lake Erie for August of 2000-2015. The total mass of phosphorus was conserved over the simulation period. In addition to photosynthetically active radiation (PAR), the biophysical model was forced by temperature and vertical turbulent diffusivity from the Lake Erie Operational Forecasting System model (LEOFS, https://tidesandcurrents.noaa.gov/ ofs/dev/leofs/leofs.html), an application of the Finite Volume Community Ocean Model (Chen et al. 2003;Anderson et al. 2015). Biophysical model simulations differed by station and year, using water column depth, temperature, turbulent diffusivity, and PAR from LEOFS model output specific to each station and year. Meteorological forcing for the hydrodynamic model (wind velocity, air temperature, humidity, solar radiation) were obtained by interpolation from meteorological stations surrounding Lake Erie, as described by Beletsky et al. (2003). Phytoplankton and zooplankton were initialized at low concentrations and allowed to approach steadystate with ambient light, temperature, and nutrient availability over the month of January. Mussel biomass was initialized at a level that resulted in little net change in mussel biomass over the 12-month simulation, based on field observations of relatively stable mussel biomass at deep stations in eastern Lake Erie (EPA GLNPO 2007-2015 Dreissena data). For our analysis of the influence of food availability on mussel growth, we assumed that food availability to the mussels Ratio of phytoplankton biomass predicted by biophysical model including vs. excluding the presence of mussels (d). Grazing reduces phytoplankton throughout the water column during mixing periods and in the hypolimnion but enhances epilimnion phytoplankton by indirectly increasing nutrient levels remaining after the spring mixing period via a top-down cascade. Mussel size data from Dermott and Munawar (1993), Mills et al. (1993), Roe and MacIsaac (1997) was proportional to the phytoplankton model variable and rescaled these depth-specific food levels to fall between 0 (minimum) and 1 (maximum). Finally, to illustrate the effect of mussel grazing on intraannual phytoplankton dynamics, we examined the ratio of phytoplankton biomass in model simulations of our offshore sites with vs. without observed levels of mussel biomass in eastern Lake Erie.

Environmental conditions during experiments
At the beginning of the 2015 experiment a weak thermal stratification developed which later strengthened throughout the experiment until the end of October when homothermy became established (Supporting Information Fig. S2). Temperatures at the littoral sites were consistently higher throughout the season than in the profundal sites (average for the season for three sites: 7.50 6 0.178C in profundal vs. 16.00 6 0.358C in littoral sites, mean 6 standard error here and throughout unless otherwise noted; Table 1) and gradually increased to a maximum of 258C in August-September. However, there were several occasions during seiche events when the temperature at the littoral sites rapidly dropped and was almost as low as the temperatures in the profundal sites (Supporting Information Fig. S2). Within treatments, temperatures were highly consistent throughout the season.
Although our 2016 study started at nearly the same time as in 2015, stratification was established 1 month later (Supporting Information Fig. S2). On the bottom and 0.7 m above the bottom temperature was low and gradually increased from 58C to 78C throughout the season. Cages suspended 25 m above the bottom (26-30 m below the surface) were largely located within the thermocline and experienced the largest temperature variations. Cages suspended 35 m above the bottom (16-20 m below the surface) were located above the thermocline, where temperature gradually increased to a maximum exceeding 258C in August and then declined. As in 2015, temperatures were consistent within treatments at similar depths, and there were several occasions when the temperature in the upper-most cages rapidly dropped almost as low as in the profundal due to seiche events (Supporting Information Fig. S1). In addition to decreasing temperatures, seiche activity also temporarily exposed these cages to water with the lower phytoplankton concentrations characteristic of the hypolimnion.
Simulations of phytoplankton distribution in our biophysical model were generally consistent with in situ measurements and illustrate the ephemeral nature of phytoplankton abundance in profundal habitats. In particular, phytoplankton distribution across depths qualitatively matches in situ fluorescence profiles at profundal stations (Fig. 5c). As highlighted by Rowe et al. (2017) in Lake Michigan, mussel grazing reduces phytoplankton densities in the hypolimnion, especially immediately following the onset of thermal stratification (Fig. 4d). Additionally, mussels also reduce food availability to the profundal benthos during key spring and fall periods of deep mixing. By contrast, littoral areas experience a relatively consistent abundance of food throughout the growing season (Fig. 5a) because the surface mixed layer and euphotic zone are connected to the bottom by vertical turbulent mixing. Here, mussels consistently reduce resources throughout the growing season, but do not deplete resource availability in benthic littoral habitats to very low levels. We also note that grazing enhances simulated epilimnion phytoplankton biomass through an indirect effect: because mussels consume phytoplankton during spring mixing periods, more unused nutrients remain in the water column as stratification sets in, promoting phytoplankton growth. These contrasting effects of mussels on hypolimnetic vs. epilimnetic phytoplankton are consistent with observed, diverging trends in mid-summer turbidity in near-bottom vs. near-surface samples (Fig. 4a). However, we emphasize that an increase in epilimnetic phytoplankton is an intra-annual effect of mussel grazing: our comparison of mussel effects does not account for the inter-annual decline of phytoplankton biomass commonly observed in waterbodies after Dreissena invasion (reviewed in Karatayev et al. 1997Karatayev et al. , 2015Higgins and Vander Zanden 2010).

2015: Mussel demography in littoral vs. profundal
Overall, mussel growth was size-dependent and was much faster in the littoral compared to profundal habitats (3.91 6 0.26 mm in littoral sites 4-6 vs. 2.15 6 0.15 mm in profundal sites 1-2; Table 1). While littoral growth rate and temperatures were highly consistent among sites, in profundal habitats growth was higher at the shallowest site 3 compared to the other sites. Conditions at this site were less representative of the other profundal sites due to the fact that this site began experiencing warmer temperatures in early October, approximately a month before other sites (Supporting Information Fig. S1). Consequently, faster growth may be explained by stronger effects of the fall bloom, consistent with our biophysical model results suggesting that most food in the profundal habitats is available to mussels during only brief mixing periods (Fig. 5a,b; see below), or by the contact of the deep chlorophyll layer with the lake bottom, which can promote mussel growth (Malkin et al. 2012). In either case, these conditions do not appear to be representative of other profundal habitats below 50 m depths. Dreissena mortality in cages over the whole period of experiment was low, but somewhat higher in littoral (3.9%, 8 died out of 281) vs. profundal habitats (0.4%, 1 died out of 285). Dead individuals were excluded from the growth rate estimates.
In 2015, we also noticed strong recruitment in littoral habitats. All parts of our mooring system, including cages, lines, cinderblocks, and even mussels inside cages were covered by the young-of-the-year Dreissena spp. (Fig. 2). The density of settled mussels obtained from the line varied among sites from 16,987 m 22 to 118,800 m 22 . Zebra mussels (D. polymorpha) dominated young-of-the-year newly settled mussels, consistently forming from 64% to 72% of all dreissenids. However, newly settled quagga mussels had higher maximum length and were significantly larger than zebra mussels. In contrast, in profundal habitats in 2015 we did not find a single young-of-the-year Dreissena on any parts of our mooring system, including cages, lines, and cinderblocks located on or 0.7 m above the bottom.

2016: Drivers of mussel demography
Dreissena recruitment was strongly depth dependent and was likely limited by very low food densities in the hypolimnion. As in 2015, we found no settlement on or near the bottom, while all parts of the mooring system, including, cages, lines, and floats suspended in the water column (25 m above the bottom and higher) were completely covered with young-of-the-year Dreissena in all three sites (Figs. 2, 6c). Large numbers of recently settled mussels were even found inside each of the three upper cages (1250-2870 mussels per cage). On floats, lines and cages closer to the surface (16-20 m below the surface) mussels density often exceeded 100,000 m 22 and declined with depth. However, quagga mussels were dominant (> 80%) in all 2016 recruitment samples. Forward stepwise Ridge regression selected relative food as the strongest predictor of recruit densities (multiple R 5 0.93, b 5 0.58, p 5 0.007), indicating that recruitment was more closely related to relative food availability than temperature. Accordingly, a sharp decline in recruitment below 25 m depths was also broadly consistent with patterns in relative food availability predicted by the biophysical model throughout the growing season (Fig. 6a), and in particular with the very low phytoplankton abundance under stratified conditions in the middle of the growing season (when most settlement is likely to occur, reviewed in Karatayev et al. 2014a;Fig. 5b).
Similarly to recruitment, mussel growth was much higher at shallower depths, and was affected by temperature (degree days), relative food availability, and the depletion of food within a near-bottom boundary layer. Compared to 2015 littoral habitats, mussel growth was almost twice as fast in the 2016 warm-water treatments, which coincided with the deep chlorophyll layer (Fig. 5b), and consequently experienced greater phytoplankton densities. In addition, they were not exposed to any boundary layer effects because cages were suspended in the water column. At the same time, growth in profundal treatments was generally lower because the spring bloom occurred earlier in 2016 vs. 2015 (with food in profundal habitats rapidly declining in day of year 100 vs. 130; Fig. 5a vs. 5b), and thus was largely missed in our 2016 profundal treatments. Bayesian regression analysis found relative food availability from the biophysical model to strongly predict von Bertalanffy growth rates (K) of mussels throughout the water column (mean of the posterior distribution: 0.49, 2.5 th and 97.5 th percentiles: 0.27 and 0.71, respectively), although the effect of degree days was also significant and appeared similar in magnitude (mean: 0.55, 2.5 th and 97.5 th percentiles: 0.43 and 0.67). Additionally, we found evidence of reduced mussel growth within a 0.7 m layer near the bottom of the lake at the three profundal sites in 2016. In particular, mean mussel growth was slowest in cages located directly on the bottom, not exceeding 1 mm for the season across all sites (Table 1) and nearly eightfold lower than at the shallowest (16-20 m) depths (0.97 6 0.01 vs. 7.46 6 0.66 mm). By contrast, growth increased by 44% in cages located just 0.7 m above the bottom compared to cages directly attached to the cinderblocks. Consistent with this pattern, Bayesian regression identified the presence of cages on the bottom as a strong negative effect on growth rates (mean: 20.48, 2.5 th and 97.5 th percentiles: 20.85 and 20.11). Again, Dreissena mortality during the experiment was very low (only two mussels died across all cages).

Long-term dynamics in mussel food resources
Analysis of the U.S. EPA GLNPO long-term data from the stations > 40 m depths in the eastern basin of Lake Erie revealed that the seston concentration dynamics from 1996 to 2006 was similar in the surface and near bottom layers (Fig. 4a,b). Over this period, average surface and bottom seston during the spring were not meaningfully different (< 9%; 0.78 6 0.20 mg L 21 vs. 0.85 6 0.22 mg L 21 near surface and bottom, respectively; Student's t 5 23.11, t-test). Likewise, patterns of summer seston concentration dynamics in 1996-2006 near the surface and bottom were similar, although overall differed more strongly compared to spring (37%; 0.40 6 0.03 mg L 21 vs. 0.29 6 0.03 mg L 21 near surface and bottom, respectively; t 5 4.08). After 2006 we found a threefold difference between surface and bottom summer seston concentrations ( Fig. 4b; 0.59 6 0.05 mg L 21 vs. 0.20 6 0.01 mg L 21 near surface and bottom, respectively; t 5 9.60; mean across 2007-2013), suggesting that depletion of food resources in the hypolimnion during stratification is substantially increased in the presence of high mussel biomass. The decline in near bottom summer seston concentration coincides with a strong increase in average mussel length in profundal habitats, compared to the earlier in the invasion (8.26 6 1.21 mm in 1992-1998 vs. 23.01 6 0.54 mm in 2008-2015, t 5 211.1, p < 0.001, Fig. 4c) due to the lack of young (< 5 mm) mussels. By contrast, mean mussel length in littoral habitats did not change over this time period (7.80 6 0.63 mm vs. 12.05 6 2.78 mm in 1992-1998 and 2008-2015, respectively; t 5 21.49, p 5 0.27; Fig. 4), reflecting a continuous presence of small recruits. Additionally, mean mussel length did not differ among sites in littoral and profundal habitats earlier in the invasion (1992)(1993)(1994)(1995)(1996)(1997)(1998), when an abundance of recruits was found in both habitats (7.80 6 0.63 mm vs. 8.26 6 1.21 mm in littoral and profundal sites, respectively; t 5 20.34, p 5 0.75). Indeed, since 2009 we did not record successful recruitment and growth of young-of-theyear Dreissena at U.S. EPA long-term monitoring stations in the eastern basin located at depths over 40 m (Fig. 1). Those years when the settlement did occur, recruits did not seem to survive to the following year.  Table 1). Note the log scale in (c).

Discussion
We found a strong gradient of declining dreissenid recruitment and growth with increasing depth, following declines in food availability and temperature. Low growth, large average size, and the lack of recruitment imply much higher longevity and much lower biomass turnover of mussels in deep vs. shallow areas. Consequently, the bulk of dreissenid biomass found in profundal habitats of deep lakes or lake basins (e.g., lakes Michigan, Ontario, and Huron) may have very limited ecological impacts on epilimnetic resources, compared to shallower warmer areas where mussels have high growth and production and therefore higher ecosystem impacts through their feeding and filtering activities. Our analysis of long-term trends in seston concentrations (turbidity) 1 m above the bottom and of our biophysical modeling results indicate that this effect is driven by depletion of a limited food supply in profundal habitats during long periods of summer stratification when phytoplankton production in the euphotic zone is not efficiently transported to the bottom by vertical mixing. At the same time, slower growth and higher longevity also explain the slower rates of biomass increase over time exhibited by profundal dreissenid populations in the Great Lakes (e.g., Rowe et al. 2015b). Finally, limited recruitment in the profundal habitats owing to competition over limited food resources may ultimately limit mussel abundance and, if sufficiently strong, may produce long-term fluctuations in quagga mussel populations and their ecosystem impacts.

Drivers of demography at depth
In littoral habitats, mussels grew nearly twice as fast as in most profundal habitats, an effect which our study likely underestimated by omitting the effects of food depletion within a near-bottom boundary layer because all cages in 2015 were suspended about 0.7 m above the bottom. Suspension of cages in the water column above the benthic layer resulted in elevated growth compared to animals on the bottom, where food resources are depleted due to mussel filtration, a result found for both fresh and marine mussels (Wildish and Kristmanson 1984;Fr echette and Bourget 1985;Peterson and Beal 1989;Yu and Culver 1999;Karatayev et al. 2006;Malkin et al. 2012). In the Great Lakes specifically, seston-depleted benthic boundary layers extending 1-2 m above dreissenid mussel beds have been reported previously (Ackerman et al. 2001;Liao et al. 2009); these layers were maintained by mussel consumption of particles and a downward turbulent flux of particles into the benthic boundary layer from waters above. However, the extent to which this phenomenon affects mussel populations has not been examined. In 2016, we found that at every site mussels suspended at 70 cm above the bottom grew faster than those on the bottom. Consequently, our 2015 estimates of the extent to which growth rates are reduced in profundal habitats compared to shallow areas appear to be conservative.
Gradients in temperature, food, and turbulent mixing among littoral and profundal habitats are most likely responsible for observed differences in mussel growth. Many studies have observed increased growth rates at higher temperatures in quagga (MacIsaac 1994;Karatayev et al. 2011b) and especially in zebra mussels (reviewed in Karatayev et al. 2006Karatayev et al. , 2011b; similarly, the importance of food as a limiting factor for Dreissena spp. growth is well known (Dorgelo 1993;Sprung 1995;Jantz and Neumann 1998;Schneider et al. 1998;Horvath and Lamberti 1999). Vanderploeg et al. (2010) found that quagga mussel clearance rate decreased linearly with temperature below 78C. In addition to the effect of temperature on biological process rates, temperature is associated with the vertical turbulent diffusivity because the thermocline suppresses turbulence; thus, when the surface mixed layer is in contact with the bottom, greater food availability and higher temperature may occur simultaneously (Edwards et al. 2005;Liao et al. 2009;Waples et al. 2017).
Because food concentrations often differ between zones with different temperatures and circulation patterns, it can be difficult to tease apart the relative importance of each factor. Indeed, studies suggesting that dreissenid growth does not occur at low temperatures (Alimov 1974;Mackie 1991;Jantz and Neumann 1992;Morton 1969a,b) have been based on field observations during the winter, when low temperatures are coincident with low phytoplankton abundance. Our Bayesian regression analysis overcomes this limitation by synthesizing results from both lab (Karatayev et al. 2011b) and field (this study) experiments, and highlights that both factors regulate mussel growth. Although our results suggest that the effect of food on growth are similar in magnitude to that of temperature, this estimate may be conservative because, near-bottom growth was almost twice as fast in 2015 as in 2016, when our profundal treatments spanned the spring bloom to a greater degree (Fig. 5a vs. 5b) but mean temperatures were very similar (Table 1). Thus, future lab studies on the effect of food levels on mussel growth would greatly aid in better resolving the independent and interactive effects of food and temperature on mussel growth.
In contrast to growth, we find that mussel recruitment was primarily determined by food availability. In particular, quagga mussel recruitment occurred both above and below the mid-summer thermocline, but we did not find a single recently settled Dreissena on our mooring system below 30 m during 2 yr of our experiments (e.g., Fig. 6c). Similarly, the near absence of dreissenid recruitment at any of four long-term monitoring EPA profundal benthic stations in the eastern basin across 7 yr of observations ( Fig. 1) generally mirrors observed patterns in near-bottom turbidity, a proxy for mussel food availability, in profundal habitats (Fig. 4b).
Recruitment must have occurred early in the invasion to produce the existing population of large mussels, but subsequently ceased when competition for limited food resources Karatayev et al. Dreissena demography in stratified lakes became too great, even though there was likely little change in temperature of the profundal benthos over the same period. Based on our field experiment and long-term population demography, we therefore hypothesize that successful Dreissena recruitment is limited by food availability rather than temperature or the presence of larvae and may be inhibited by competition of juveniles with large mussels. Our analysis also indicated that temperature and food gradients between littoral and profundal habitats are in turn both driven by summer thermal stratification and concomitant food depletion in the profundal zone. In general, survival and growth of quagga mussels as benthic filter feeders is dependent on the particle flux from the water column to the bottom. Dreissenids are known to reduce food resources (e.g., seston and chlorophyll concentrations) at both lakewide scales (reviewed in Karatayev et al. 1997Karatayev et al. , 2007Higgins and Vander Zanden 2010;Pothoven and Fahnenstiel 2013;Rowe et al. 2015b) and more locally within lake basins during thermal stratification (MacIsaac et al. 1992;Ackerman et al. 2001;Edwards et al. 2005). For example, Rowe et al. (2015b) showed reduced satellite derived chlorophyll a (Chl a) concentration in Lake Michigan during periods of deep mixing, spatially associated with areas of high mussel biomass around the 30-m to 50-m depth contours. Pothoven and Fahnenstiel (2013) found that Chl a in the hypolimnion (bottom 20 m) of Lake Michigan declined 63% and 54% between the pre-and post-quagga mussel invasion periods in early and late summer, respectively. Our biophysical model reinforces these results, indicating reduced food abundance due to filter feeding at deep stations during periods of deep mixing when food availability is greatest, and in the summer hypolimnion when food is very limited (Fig. 5a; Rowe et al. 2015aRowe et al. , 2017. Moreover, given that long-term declines in profundal seston concentrations in eastern Lake Erie have become apparent only in recent years, we suggest that intense food depletion in the profundal zone is associated with the recent climax in profundal Dreissena biomass (Fig. 4c).

Implications for ecosystem impacts of dreissena
Being extremely important players in Great Lakes ecosystems, dreissenids greatly enhance benthic-pelagic coupling, increasing water clarity, decreasing phytoplankton and primary production in the water column (reviewed in Fahnenstiel et al. 2010;Pothoven and Fahnenstiel 2013;Karatayev et al. 2015), which can trigger a suite of connected changes in the entire ecosystem, including changes in benthic communities and fish (Hoyle et al. 1999;Lozano et al. 2001;Pothoven et al. 2001;Madenjian et al. 2002;Nalepa et al. 2009a,b). However, the ecological impact of dreissenids is reflected in their growth and productivity rather than density alone. Mussels in profundal habitats may grow at least twice slower compared to littoral sites, and individuals at 50 m depth therefore reach the average adult size at an older age. Given that mean lengths of multiyear adult mussels in profundal habitats are similar or larger than those in littoral habitats (Fig. 1), this difference in growth implies much higher survival rates and longevities of mussels in deeper areas (potentially due to lower metabolic costs, physical wave stress, and predation). Slow growth and high longevity, in turn, translate into much lower biomass turnover (i.e., productivity) for mussels in profundal habitats.
Habitat differences in mussel productivity implied by our results have key implications for other Great Lakes (e.g., Michigan, Ontario, and Huron), where both overall dreissenid biomass and its recent increase have occurred primarily in the extensive profundal habitats (Rowe et al. 2015b;Nalepa and Baldridge 2016), as well as for lakes with long periods of summer stratification more broadly. First, observed declines in mussel growth rates with depth and in profundal vs. littoral benthic habitats reflect a strong gradient in overall dreissenid impacts on pelagic resources because mussels in profundal habitats have access to (and consequently impact) epilimnetic primary production only during brief mixing periods. This result markedly contrasts predictions of mussel impacts based on dreissenid biomass alone, which is nearly twofold higher in profundal vs. littoral habitats in our system (Karatayev et al. 2014b). Additionally, our study provides experimental verification to the biophysical modeling results of Rowe et al. (2017) predicting that mussels in profundal habitats of Lake Michigan are strongly resource-limited. A second implication of our experimental and modeling results is a strong effect of mussel grazing on food availability to organisms in profundal habitats. In particular, mussels reduce near-bottom food levels not only during stratified periods (possibly by > 50%; Fig. 4d; Pothoven and Fahnenstiel 2013;Rowe et al. 2017), but also during spring and fall mixing periods, which represent key resource pulses to benthic organisms in the profundal zone.
Increased mussel longevity (> 7 yr) implied by the strong difference in growth rates in littoral and profundal habitats also explains the slower dynamics of Dreissena densities in profundal vs. littoral habitats observed throughout the Great Lakes. For example, quagga mussels in the profundal habitats of eastern basin of Lake Erie experience more stable longterm density dynamics compared to the much shallower western basin (33 lower coefficient of variation in densities; 2003-2013 U.S. EPA GLNPO data). More broadly, initial peaks in zebra and quagga mussel populations both in North America and Europe occur several years after their invasion, and are generally followed by a limited decline (reviewed in Karatayev et al. 2015). Although similar dynamics have been observed in the nearshore zones (0-30 m) of lakes Erie, Michigan, and Ontario, offshore densities and biomass stabilized only in Lake Erie, which has the smallest profundal zone and was invaded by quagga mussels earlier than other Great Lakes (Karatayev et al. 2014b). In deeper (> 90 m) profundal habitats of lakes Michigan and Ontario quagga mussels density continue to slowly increase Rowe et al., 2015a;Nalepa and Baldridge 2016), suggesting that mussels have not yet exhausted their food resources and reached their carrying capacity in the deepest areas of these lakes.
Finally, our results indicate the presence of density dependence of mussel recruitment in the profundal areas of eastern Lake Erie, which may lead to long-term fluctuations in dreissenid populations. In general, strong density dependence may induce generational fluctuations in populations, as abundant cohorts inhibit recruitment in future years. The spatial extent of such fluctuations in our system will depend on the scale at which mussels compete for food resources in profundal habitats. For example, if high mussel densities only deplete seston (and therefore recruitment) at small scales, population fluctuations will be relatively local and average out over larger (i.e., lakewide) scales. By contrast, if food depletion occurs synchronously throughout the whole hypolimnion, successful recruitment in the eastern basin of Lake Erie may be possible only after the bulk of old mussels dies and more food resources become available for juveniles.

Conclusions
We find a drastic gradient in quagga mussel demography, including recruitment, growth and likely longevity between littoral and profundal benthic habitats which results from prolonged summer thermal stratification, a distinct property structuring the environment of numerous temperate lakes. These results imply reduced Dreissena productivity in profundal habitats and reflect much weaker, periodic grazing impacts on epilimnetic primary production, consistent with predictions of recent ecosystem-scale biophysical models. Thus, our results explain how the ecosystem impacts of one of the most aggressive freshwater invaders are likely to vary among lakes depending on morphometry (i.e., size of profundal areas) and the strength and duration of summer stratification. More broadly, our study highlights the degree to which the ecological impacts of invasive species may vary both across ecosystems and across habitats within the same environment. Greater emphasis on identifying and characterizing such spatial heterogeneity is key to understanding and managing the effects of species invasions.