Introduction

Climate change, land use change and species invasions are among the most impactful stressors that threaten the Laurentian Great Lakes (hereafter referred to as “Great Lakes”) (Allan et al. 2013). Recent climate trends show that the Great Lakes surface water temperatures are warming faster than the average global rate (0.34 °C decade−1), with Lake Michigan’s surface temperature increasing at a rate of 0.2–0.42 °C decade−1 (O’Reilly et al. 2015; Collingsworth et al. 2017). Empirical studies and models indicate that climate warming in deep temperate lakes has shortened durations of ice cover and isothermal periods, while advancing the onset and lengthening the period of summer stratification (Brooks and Zastrow 2002; Winder and Schindler 2004; Shimoda et al. 2011; Mason et al. 2016). Climate change may also affect the biological character of the Great Lakes by altering the timing of phytoplankton blooms (Winder and Schindler 2004), affecting the availability and transport of nutrients throughout the water column in deep lakes (Brooks and Zastrow 2002), and increasing the volume of favorable thermal habitat available to most fish species thereby affecting growth and consumption rates (Collingsworth et al. 2017). Meanwhile, more frequent and intense runoff events due to changing precipitation patterns (Michalak et al. 2013; Bosch et al. 2014), paired with recent trends in urbanization and population growth (Pijanowski and Robinson 2011), threaten to increase non-point and point-source nutrient inputs in areas of the Great Lakes already struggling with eutrophication, such as Green Bay, Saginaw Bay, and Lake Erie’s Western Basin (Stow et al. 2015; Choquette et al. 2019).

The potential effects of climate change and land use change on the Great Lakes’ ecosystem should be understood in the context of historical changes in which species invasions have been a primary factor (Mills et al. 1993; Madenjian et al. 2008; Cuhel and Aguilar 2013). Among the 180+ non-indigenous species in the Great Lakes, zebra mussels Dreissena polymorpha and quagga mussels D. rostriformis bugensis (hereafter collectively referred to as ‘dreissenid mussels’) have caused some of the most substantial ecosystem changes (Fahnenstiel et al. 2010; Vanderploeg et al. 2010). In conjunction with nutrient load reductions, dreissenid mussels have contributed to declines in phytoplankton biomass and affected spatial and temporal patterns of productivity in Lake Michigan (Fahnenstiel et al. 2010; Mida et al. 2010; Vanderploeg et al. 2010; Rowe et al. 2017). As managers work to prevent future species invasions, there is a need to understand how the Great Lakes’ susceptibility to future invasions has changed and will change over time in response to the individual and collective effects of these stressors.

Anthropogenic stressors and biological invasions may alter biotic and abiotic properties of aquatic ecosystems that previously acted as ecological barriers to the establishment of invasive species (Hellmann et al. 2008; Rahel and Olden 2008). However, the multitude of direct and indirect effects of stressors, such as climate change, make it difficult to predict how invasion risk will change (Rahel and Olden 2008). As a further complication, multiple stressors can interact in sequential, additive, synergistic, or antagonistic ways (Kelly et al. 2017; Sarà et al. 2018; Smith et al. 2019). The uncertain effects of multiple stressors on ecosystem ‘invasibility’ motivates risk assessment research that assesses future invasion vulnerability based on current trajectories of multiple anthropogenic stressors (Smith et al. 2012).

Bighead Carp Hypophthalmichthys nobilis and Silver Carp Hypophthalmichthys molitrix (collectively ‘bigheaded carps’ [BHC]) are two species of high concern due to their proximity to the Great Lakes and the potential impact they could have on the ecosystem. These highly invasive, planktivorous cyprinids have established prolific densities in parts of the Mississippi River basin, including the Illinois River (Garvey et al. 2015), which is considered the primary source of invasion into the Great Lakes due to its hydrological connection with Lake Michigan via the Chicago Area Waterway System (Kolar et al. 2007). In North America, the intensive grazing of BHC on plankton has altered size structure and composition of zooplankton communities (Sass et al. 2014) and reduced algal biomass (Tumolo and Flinn 2017), which consequently has contributed to the declining growth condition and abundance of several resident planktivores (Irons et al. 2007; Pendleton et al. 2017). Because BHC feed primarily on plankton, they are capable of affecting the growth rates of most fish species, not just planktivores, during their plankton-dependent larval stage (Fletcher et al. 2019). The introduction of these efficient filter-feeders into Lake Michigan could put more stress on prey fish populations, which have been declining over the past few decades (Madenjian et al. 2012). However, the realization and severity of BHC impacts depends on their ability to survive and grow in a colder, less productive environment than they typically encounter across their global range.

Earlier skepticism regarding the plausibility of BHC to be able to survive and grow in the Great Lakes has been addressed through spatially-explicit models of BHC growth rate potential (GRP) (Cooke and Hill 2010; Anderson et al. 2015, 2017; Alsip et al. 2019). In our previous study, Alsip et al. (2019), we built on the work of Cooke & Hill (2010) and Anderson et al. (2015, 2017) by demonstrating how bioenergetically-suitable habitat for BHC is not only present in Lake Michigan, but more expansive than originally thought if the fish leveraged their characteristically-flexible diet and utilized subsurface habitat. High quality habitat for BHC growth was concentrated in a few highly productive areas, including Green Bay, whereas offshore areas were only capable of supporting weight maintenance for Bighead Carp. While Lake Michigan’s cold and food-scarce waters are not an insurmountable barrier for BHC growth, there is a need to understand how alterations to the lake’s productivity and thermal regime as a result of anthropogenic activity affects the lake’s suitability for BHC.

The purpose of this paper is to investigate the mechanisms by which climate, nutrient loads, and invasive dreissenid mussels individually and interactively affect Lake Michigan’s bioenergetic suitability for BHC. We build on our previous work in Alsip et al. (2019) by applying a GRP model (Brandt et al. 1992) to quantify habitat quality using scenario datasets from a 3-dimensional biophysical model of Lake Michigan developed by Rowe et al. (2017). GRP provides a quantitative measure of fish habitat quality by translating measurements of environmental conditions (prey biomass and water temperature) into fish growth (g g−1d−1) based on bioenergetics. Rowe et al. (2017) generated the scenario datasets with their biophysical model to investigate the influence of meteorology, tributary phosphorous loads, and dreissenid mussel filtration on primary production. The availability of these simulated environmental data provided an opportunity to evaluate the effect of these stressors on BHC habitat suitability. For individual stressor effects, we hypothesized that the negative impact of mussels on suitable BHC habitat would be spatially and temporally discrete—being confined to isothermal periods and to areas where mussels are in contact with the surface mixed layer (SML). We also hypothesized that a warmer climate and increased nutrient loads would limit the negative effects of mussel filtration on BHC habitat suitability. Specifically, we expected a warmer climate and increased nutrient loads to increase the duration of the growing season, the quantity of suitable habitat, and the annual growth potential of BHC.

Methods

Study site

Lake Michigan is a meso-oligotrophic lake at temperate latitudes (Fig. 1). It has a surface area of 57,800 km2, a mean depth of 85 m, a maximum depth of 282 m, and average summer surface temperatures that reach 21–22 °C (NOAA Great Lakes CoastWatch Program https://coastwatch.glerl.noaa.gov/). Lake Michigan’s biotic and abiotic environment is spatially heterogeneous and dynamic (Rowe et al. 2017). The lake is dimictic—mixing in the spring and fall and thermally stratifying in the summer and winter. Annual maximum ice cover ranges from 12.4 to 93.1% (1973–2019 mean = 40%) (https://www.glerl.noaa.gov/data/ice/#historical).The formation of a thermocline during summer stratification divides the water column into three ecologically distinct zones: an epilimnion, metalimnion, and hypolimnion. Deep chlorophyll layers occur during summer at an average depth of 30 m (Bramburger and Reavie 2016).

Fig. 1
figure 1

Map of Lake Michigan (a), showing the spatial domain of Finite Volume Community Ocean Model – General Ecosystem Module (FVCOM-GEM; white area), bathymetry (50-m contours), bordering states (bolded names), Chicago (high potential introduction location; circle), inflows of the 38 tributaries (triangles) included in FVCOM-GEM with the largest phosphorus contributors labeled by name, and the mid-depth location near Muskegon, MI (square). Enlarged area of southeastern Lake Michigan (b), showing the hydrodynamic model grid, the mid-depth location near Muskegon (square labeled M45), and the location of tributary mouths (triangles)

Full size image

Growth rate potential model

We used a spatially explicit model of fish GRP to evaluate how climate, nutrient loads and dreissenid mussels affect habitat suitability for BHC. GRP model development and structure is fully described in Alsip et al. (2019). The model integrates outputs from a three-dimensional biophysical model, a foraging model to estimate consumption rate (g g−1d−1), and a bioenergetics model to estimate GRP (g g−1d−1) for a non-swimming and non-reproducing fish. The GRP model is constrained by species-specific physiological parameters and is driven by habitat conditions (i.e. temperature and prey concentrations). Inputs to the foraging model and bioenergetics equations in the GRP model were output from the spatially explicit biophysical model. Prey inputs included phytoplankton, zooplankton, and detritus, the latter of which was characterized with the energetic quality of dreissenid biodeposits as determined by Anderson et al. (2016) and adjusted by an assimilation coefficient derived from that same study. We used prey energy density values of 2600 J g−1 wet mass, 2512 J g−1 wet mass, and 127.3 J g−1 wet mass for phytoplankton, zooplankton, and detritus, respectively (Anderson et al. 2015, 2016, 2017). All model parameters and equations can be found in the supplement to this manuscript (Tables S1.1–S1.3) and additional information on model mechanics and parameterization can be found in Alsip et al. (2019) supplementary information. All simulations were coded and run in R version 3.5.1 (https://CRAN.R-project.org). Each component of the GRP model is discussed briefly below.

Bioenergetics model

We used the Wisconsin Fish Bioenergetics 3.0 model (Hanson et al. 1997), which uses a mass balance approach that estimates growth rate (G, g g−1d−1) of an individual fish by subtracting respiration (R), egestion (F), excretion (U), and specific dynamic action (S) from estimates of consumption (C):

$$ G = C - \left( {R + F + U + S} \right) $$
(1)

Foraging model

We calculated C by taking the minimum value of two consumption estimates: maximum consumption based on mass and temperature (Cmax, Supplementary information, Table S1.2) and foraging-based consumption (\(C_{FR} )\). Cmax is determined by the bioenergetics equation for consumption whereas \(C_{FR}\) is a function of temperature (\(f\left( T \right))\), prey concentration (g L−1), and filtration rate (FR; L d−1), which itself is a function of fish mass W (g) and foraging hours (t) (from Smith, 1989):

$$ C_{FR} = \left( {FR *\frac{{\left( {PP + ZP + Det.} \right)}}{W}} \right) * f\left( T \right) $$
(2)
$$ FR = 1.54 * W^{.713} * t $$
(3)

We set t = 24 hours for each simulation. We then multiplied the minimum value between \(C_{FR}\) and Cmax by a prey-to-predator energy density (ED) ratio to calculate C (g g−1d−1):

$$ C = min\left( {Cmax, C_{FR} } \right)*\frac{{ED_{Prey} }}{{ED_{Carp} }} $$
(4)

Biophysical model

The three-dimensional, heterogeneous environment was defined by prey concentrations and water temperatures simulated by the Lake Michigan Finite Volume Community Ocean Model–General Ecological Module (FVCOM-GEM, Fig. 1) (Rowe et al. 2015, 2017). FVCOM is a 3-dimensional, hydrodynamic numerical model that predicts currents, temperature, and water levels driven by external physical forcing including surface wind stress, and heat flux (Chen et al. 2006). The unstructured grid and terrain-following sigma vertical coordinate of the model allows for accurate representation of complex coastline morphology. FVCOM includes a General Ecological Module (GEM), which allows for flexible representation of the lower food web (Ji et al. 2008). FVCOM was applied to Lake Michigan using 20 sigma layers of uniform thickness, and an unstructured grid consisted of 5,795 nodes and 10,678 elements, with element side lengths of 0.6 to 2.6 km near the coast and 4.5 to 6.8 km near the center of the lake (median 3.1 km) (Rowe et al., 2015). Rowe et al. (2017) implemented GEM as a phosphorus-limited nutrient-phytoplankton-zooplankton-detritus (NPZD) model that simulated lower food web biomass and productivity and included a dreissenid mussel (benthic filter feeder) compartment. Phosphorus loads from 38 tributaries were included in FVCOM-GEM. The geographic scope of our GRP model was confined by the boundary of FVCOM’s spatial grid, which included Lake Michigan and Green Bay, but not upstream tributaries or drowned river mouths (Fig. 1) (Rowe et al. 2015, 2017). Model development and skill assessment was reported by Rowe et al. (2015, 2017) with an additional assessment for Green Bay and Muskegon reported by Alsip et al. (2019).

Model scenarios

We ran our GRP model with the nine biophysical model scenarios described by Rowe et al. (2017) (Table 1). Scenarios were defined by the combination of three attributes: presence/absence of dreissenid mussels, annual tributary phosphorus load levels, and meteorological forcing. Tributary phosphorous loads varied between 0 metric tons per annum (MTA) (“No Loads”), reference loads that represented 2010 levels (“Loads” = 3300 MTA), and high loads (“High P” = 5600 MTA) that matched the target load for Lake Michigan established in the Great Lakes Water Quality Agreement and were similar to loads in the 1980s (Dolan and Chapra 2012; Rowe et al. 2017). In the “High P” scenario, total phosphorous and dissolved phosphorous concentrations in the lake were initialized with data from the EPA’s Great Lakes National Program Office’s spring survey in the 1980s. Dreissenid mussel biomass in scenarios with mussels present were initialized from a 2010 benthic survey (Rowe et al. 2015, 2017). We evaluated climate effects by using meteorological forcing from years with a colder-than-average winter-spring period (1997; hereafter referred to as ‘cool year’), a warmer-than-average winter-spring period (1998; hereafter referred to as ‘warm year’), and typical thermal conditions for Lake Michigan (2010; hereafter referred to as ‘reference year’) (Rowe et al. 2015, 2017). Rowe et al. (2017) generated the warm year biophysical dataset with meteorological data from 1998. Mean water temperatures in 1998 were within the range of values expected in the years 2030–2090 under a warming climate as simulated by a Coupled General Circulation Model (Boer et al. 2000; Brandt et al. 2002).

Table 1 Attributes of the biophysical model scenarios (Rowe et al. 2017) that were used as input to the growth rate potential (GRP) model in this study
Full size table

Alsip et al. (2019) ran the GRP model using simulated biophysical data representative of 2010 conditions, corresponding to the first scenario in Table 1 and is hereafter referred to as our ‘reference scenario’ representing recent lake conditions. Effects of each stressor (mussel filtration, meteorology, and nutrient loads) on habitat suitability were tested individually by comparing the reference scenario to scenarios where a given stressor was altered (Table 2). We also focused on evaluating the interactive effects of mussels with meteorology and mussels with loads to determine the capacity of mussels to influence BHC habitat suitability considering changing environmental conditions.

Table 2 Comparative analysis of scenarios to assess individual and interactive effects of meteorology (i.e., climate), mussel filtration, and nutrient loads on bigheaded carp habitat suitability
Full size table

Habitat suitability assessments

For each scenario, we evaluated the abundance, quality, duration, and distribution (hereafter collectively referred to as "habitat suitability metrics”) of suitable habitat across Lake Michigan. We defined suitable habitat as any model grid cell that could support a non-negative growth (GRP ≥ 0 g g−1 d−1, i.e. at a minimum, the carp maintains its weight), whereas habitat quality refers to the GRP value estimated for a given grid cell. We calculated the surface area and volume of suitable habitat by running the model on a monthly time step for each scenario using biophysical model data from the middle day of each month. More details on these calculations and our basic assumptions for these quantity metrics (i.e. volume and surface area) can be found in Alsip et al. (2019). Habitat duration was determined by running the model on a daily time step for the entire year and summing the total number of days with suitable habitat present in each individual grid cell. We used two different approaches to analyze habitat quality. We mapped instantaneous GRP to visualize suitable habitat dynamics and we calculated annual growth potential (or annual weight change), which was used to determine the GRP of BHC integrated over the course of one year.

Annual growth

For each scenario, we ran the GRP model on a daily time step in each cell of the model grid for a year to determine the change in weight (W) that BHC would experience in Lake Michigan. The initial mass Wi=1 was set at 4350 g and 5480 g for Silver and Bighead Carp (Anderson et al. 2015), respectively in every cell j in the biophysical model’s grid. Final mass at the end of day i in cell j defined the initial mass of the following day (Wi +1, j) and is the function of GRPij,Wij, and time step (t in days).

$$ W_{i + 1, j} = (W_{ij} *GRP_{ij} *t_{i} ) + W_{ij} $$

GRP was not calculated in cells where temperature (T) was at or below 0 °C (hereafter referred to as frozen cells) as these habitats would realistically be unavailable due to ice cover. When temperatures in frozen cells became habitable (i.e. T became greater than 0 °C), Wij in those cells was set equal to the mean Wi of non-frozen cells in the same water column or was interpolated from nearby cells via inverse distance weighting if the whole water column was frozen. For each scenario, we calculated the volume (V) weighted mean change in weight for both species in the nearshore cells and offshore cells separately at the final time step (WF):

$$ \Delta W = \frac{{\sum V_{j} *(W_{F, j} - W_{i = 1,j} )}}{{\sum V_{ j} }} $$

We defined the nearshore and offshore zones based on established bathymetric thresholds (nearshore: < 30 m; offshore: ≥ 30 m) that are commonly used to delineate biophysically-distinct habitats in the Great Lakes (Riseng et al. 2018).

Scenario-specific analyses

In addition to the previously mentioned analyses, we ran additional analyses specific to a few scenarios to capture potentially unique effects of meteorology and its interaction with mussels.

Vertical analysis: mussels and meteorology

We investigated dynamics of prey, temperature, and habitat quality in the water column at a mid-depth location (M45, max depth = 45 m) near Muskegon, MI (Fig. 1) to visualize how mussel filtration and meteorology affected the vertical and seasonal distribution of habitat. We focused on these two stressors for this analysis because meteorology drives SML depth which can affect primary productivity and controls the access mussels have to food throughout the water column (Rowe et al. 2015). For this analysis, we used scenarios 1, 2, 4, 5, 6, and 7 (Table 2).

Thermal habitat analysis

We analyzed changes to thermal habitat, habitat quality, and habitat duration within the top 0–10 m of the water column (i.e. average across cells in this depth range) in response to different meteorological scenarios. We ran the model on a daily step for the whole year and during the months of January, February, March, April, and May (hereafter referred to as winter-spring months). We focused on this depth range because this is where we expected to see the greatest change in thermal habitat and it is the range of depths at which BHC typically occupy in the Illinois River (DeGrandchamp et al. 2008; Garvey et al. 2015). We quantified thermal habitat using average lake-wide temperatures and cumulative degree days greater than 0 °C (CDD). The three biophysical model data sets we used for this analysis differed only by the meteorological forcing (cool year, reference year, and warm year). Dreissenid mussels were present and tributary nutrient loads were held constant across all scenarios (Tables 1 and 2, scenarios: 1, 4, & 7).

Results

Effect of mussels

Removing mussels from the reference scenario improved the yearly average of most habitat suitability metrics for both species except for offshore annual weight change, which was relatively unchanged from the reference scenario (Fig. 2; Table S2.1). However, annual improvements from removing mussels were small for most metrics for Bighead Carp except for nearshore final weight (+ 326 g) and habitat duration (growing days + 16 days). Between the two species, Silver Carp experienced greater relative increases in habitat quantity and duration metrics, but Bighead Carp scored much higher in all metrics than Silver Carp in all scenarios (Fig. 2). Nearshore and offshore average annual weight change remained negative for Silver Carp despite the absence of mussels. Averaging suitability metrics at the seasonal scale indicated that the impact of mussel filtration was 3.3 × (Bighead) and 1.6 × (Silver) greater on spring (March, April, May) habitat extent than it was on the annual average extent (Fig. 3). For Bighead and Silver carp, respectively, removing mussels increased the extent of suitable spring habitat by 3981.2 km2 (76.5% increase) and 562.5 km2 (77.7% increase).

Fig. 2
figure 2

Bighead Carp Hypophthalmichthys nobilis and Silver Carp H. molitrix habitat suitability metrics for each scenario evaluated. Top left: Average growing days; Midde left: Average volume of suitable habitat; Bottom left: Average extent of suitable habitat; Top right: Average annual weight change (%) from initial weight in nearshore (< 30 m); Bottom right: Average annual weight change (%) from initial weight in offshore (≥ 30 m). Initial weight was 5480 g for Bighead Carp and 4350 g for Silver Carp. Error bars indicate the range of monthly means. Reference scenario: Mussels, Loads

Full size image
Fig. 3
figure 3

Bighead Carp Hypophthalmichthys nobilis and Silver Carp H. molitrix mean seasonal suitable habitat extent in the three meteorological scenarios (cool, reference, and warm) with and without mussels. Spring: March, April, and May; Summer: June, July, and August; Fall: September, October, and November. Error bars represent seasonal minima and maxima

Full size image

Effect of meteorology on thermal regime and habitat suitability

The warm year improved CDD and habitat suitability metrics across the year when compared with the cool year (Fig. 2; Table 3), whereas differences between warm and reference meteorology were greatest at seasonal scales. Winter-spring CDD increased by 56 (Warm-Reference), 192 (Reference-Cool), and 248 degree days (Warm-Cool) and corresponded with respective increases of 10, 23, and 32 annual growing days lake-wide for Bighead Carp (Table 3); although most of the offshore areas in the southern basin increased by as much 19–36 days in the former two comparisons and 36–70 days in the latter comparison (Fig. 4). Silver Carp annual growing days increased by 2–4 days in response to warmer climates (Warm-Cool, Reference-Cool, Warm-Reference). Annual growth improved most notably in response to increased growing days and CDD in Green Bay and nearshore areas for Bighead Carp, and only in Green Bay for Silver Carp (Fig. 5).

Table 3 Thermal and suitable habitat metrics for scenarios with different meteorological forcing
Full size table
Fig. 4
figure 4

Change in total growing days in the top 10 m of the water column for Bighead Carp Hypophthalmichthys nobilis and Silver Carp H. molitrix between meteorological scenarios with mussels present (scenario 1, 4, and 7). Annual P loads were maintained at 3300 MTA for each scenario

Full size image
Fig. 5
figure 5

Change in end-of-year weights in the top 10 m of the water column for Bighead Carp Hypophthalmichthys nobilis and Silver Carp H. molitrix between meteorological scenarios with mussels present (scenario 1, 4, and 7). Annual P loads were maintained at 3300 MTA for each scenario

Full size image

Averaging across only the spring months (March–May) revealed differences in spring habitat quantity among the scenarios and highlights the period of greatest difference between warm and reference years (Figs. 2, 6, and S2.1). Relative to the corresponding reference year scenario, the warm year scenario with mussels produced an additional 4563 km2 (87.7% increase) and 84.6 km2 (11.7% increase) of suitable habitat surface area and an additional 71 km3 (107.4% increase) and 1.6 km3 (32.3% increase) of suitable habitat volume in the spring for Bighead and Silver Carp, respectively. When compared against the cool year, the warm year produced 7611 km2 (353% increase) and 527.1 km2 (187.1% increase) more suitable habitat surface area and 112 km3 (414.6% increase) and 4.8 km3 (282.8% increase) more suitable habitat volume for Bighead and Silver Carp, respectively. The early offshore expansion in May is ultimately what caused the increases in average extent and volume of spring habitat. This expansion did not occur until June for the other climate scenarios (Fig. 6).

Fig. 6
figure 6

Bighead Hypophthalmichthys nobilis and Silver carp H. molitrix growth rate potential (GRP; g g−1 d−1) during the spring–summer transition. Gray areas indicate unsuitable habitat (GRP < 0 g g−1 d−1). All scenarios in this figure have mussels present and reference year loads (3300 MTA)

Full size image

Effect of reducing Phosphorus Loads

Changing tributary nutrient loads caused substantial changes in BHC habitat suitability. The “Mussels, No loads” scenario (scenario 3) was the most restrictive for BHC growth, and only Bighead Carp achieved positive annual weight change (Fig. 2). In this scenario, Bighead Carp average final weights in the nearshore zone were reduced from the reference scenario by more than 1 kg (14.4% decrease from reference scenario, Table S2.1), and in the offshore the lack of tributary loads resulted in Bighead Carp only growing 1.86% of their initial weight over the entire year compared to 7% simulated in the reference scenario (Figs. 2, 7; Table S2.1). Interestingly, average monthly suitable habitat extent and volume over the year were only 3.6% and 7.2% lower than the reference scenario for Bighead Carp, with the greatest reduction in spring months (average reduction of 68% in extent and 70% in volume) and little change in summer or fall months (average reduction of < 1% in extent and 4–8% in volume) (Fig. 8 and S2.2). Silver Carp experienced the greatest weight loss in the nearshore and offshore in the “Mussels, No loads” scenario (Figs. 2, 9; Table S2.1). Both species experienced reduced habitat quality in the summer in this scenario, but Silver Carp habitat was completely absent in the spring and confined to Green Bay in June–October whereas Bighead Carp habitat covered most of the lake from June–November (Figures S2.3 and S2.4).

Fig. 7
figure 7

Bighead Carp Hypophthalmichthys nobilis annual growth (%) in all scenarios evaluated. Note that initial weight was 5480 g, and (c) represents the reference scenario (scenario 1, Table 1)

Full size image
Fig. 8
figure 8

Bighead Carp Hypophthalmichthys nobilis and Silver Carp (H. molitrix) mean seasonal suitable habitat extent in the three total phosphorus load scenarios (0 MTA, 3300 MTA, and 5600 MTA) with and without mussels. Spring: March, April, and May; Summer: June, July, and August; Fall: September, October, and November. Error bars represent seasonal minima and maxima

Full size image
Fig. 9
figure 9

Silver Carp Hypophthalmichthys molitrix annual growth (%) in all scenarios evaluated. Note that initial weight was 4350 g, and (c) represents the reference scenario (scenario 1, Table 1)

Full size image

Effect of increasing nutrients

The two “High P” scenarios resulted in the longest growing seasons, the greatest quantity of suitable habitat, and the highest final weights of all nine scenarios (Fig. 2). Suitable habitat extent in “High P” scenarios increased from corresponding reference load scenarios the most in spring months for both species (Average increase in spring extent: Bighead =  + 19,983 km2 and + 22,600 km2; Silver =  + 651.3 km2 and + 1255.7 km2, with and without mussels, respectively) (Fig. 8). Summer (June, July, August) and fall (September, October, November) habitat quantity for Bighead Carp did not substantially change in response to increased loads since these seasons were already saturated with suitable habitat in the reference load scenarios. However, Silver Carp habitat extent in summer and fall did increase by an average of 31,854 km2 and 32,211 km2 (summer; with and without mussels, respectively) and 34,250 km2 and 34,738 km2 (fall, with and without mussels, respectively). The increase in spring habitat quantity simulated in the “High P” scenario was due to more abundant nearshore habitat in early spring for both species and earlier expansion of offshore habitat for Bighead Carp (Figures S2.5–S2.8). Similar to habitat projected under a warm climate, suitable habitat for Bighead Carp expanded into the offshore in May in the “High P” scenarios. However, the extent of suitable Bighead Carp habitat in May surpassed that of the warm climate scenario. Silver Carp habitat expanded into the offshore in July and only did so in the “High P” scenarios (Figures S2.6 & S2.8).

Compared to the reference scenario, Bighead Carp in the “High P” scenarios grew an additional 26%–40% of their initial weight in the nearshore and an additional 24% in the offshore (Figs. 2, 7). For Silver Carp, “High P” scenarios were the only scenarios that produced substantial amounts of suitable offshore habitat (Figures S2.6 and S2.8), supported average annual growth in the nearshore (Fig. 2; grew 7.8% of initial weight without mussels and 3.3% growth with mussels), and resulted in the greatest amount of suitable habitat outside of Green Bay (Figs. 9, S2.6, and S2.8). Compared to the reference scenario, “High P” scenarios increased annual weight change for Silver Carp by an additional 10–15% in the nearshore and by 8% in the offshore. However, annual weight change remained negative in the offshore, despite the relative increase from the reference scenario and the availability of suitable habitat at certain times during the year.

Mussel impact under increased nutrient loads

While both “High P” scenarios produced the most favorable habitat conditions of all the scenarios, habitat suitability metrics responded in various magnitudes to the presence of mussel filtration. The greatest contrast between these scenarios was in habitat duration in which removing mussels increased the number of growing days by 31.4 days for Bighead Carp and 9.3 days for Silver Carp (Fig. 2; Table S2.1). The presence of mussels in the “High P” scenario reduced annual nearshore growth by 745 g for Bighead Carp and by 196 g for Silver Carp. However, differences in offshore final weights between these scenarios were negligible (scenario differences in final Bighead Carp weight < 17 g and Silver Carp weight < 3 g). Spatial and temporal patterns of mussel influence in the “High P” scenarios were similar to what was simulated in the reference scenarios in that mussels were most impactful in Green Bay and nearshore zone of the southern basin, particularly during spring and December (Figs. 8 and S2.5-S2.8).

Mussel Impacts across meteorological scenarios

Comparing between mussel and no mussel scenarios with the same meteorological year (i.e. scenarios 1 vs. 2, 4 vs. 5, and 6 vs. 7) elucidated how meteorology modulated mussel impacts on spring habitat quantity across the lake (Figs. 3 & S2.1). Mussel filtration reduced spring suitable habitat extent by 2108.1–3981.2 km2 (38.4%–97.8%) and 461.8–708.3 km2 (57.1%–251.4%) for Bighead and Silver carp, respectively. The relative magnitude of mussel-induced changes to volume was greater than that indicated by extent, but both metrics revealed that the relative negative impacts of mussels were greatest in the cool year and lowest in the warm year. Mussel impacts on annual averages of habitat quantity were less pronounced than their impact in the spring, but still notable. For the whole year, mussel filtration reduced mean suitable habitat extent and volume by an average of 1100.5 km2 (− 3% from corresponding “no mussel” scenarios) and 17.3 km3 (− 2.5%) for Bighead Carp and 423 km2 (− 22%) and 5.14 km3 (− 26%) for Silver Carp across all meteorological scenarios (Fig. 2). The average reduction in growing days across all meteorological scenarios due to mussels was 17 days (BC) and 5.8 days (SC), which corresponded with a reduction in final weights ranging between − 430.4 g to − 326.2 g for Bighead Carp and − 130.1 g to − 106.5 g for Silver Carp (Fig. 2). Differences in spatial distribution of Bighead Carp habitat quality (as measured by annual growth) were more apparent when comparing across climate scenarios than comparing across mussel scenarios (Fig. 7).

Habitat dynamics in the water column

Meteorological influence on the structure and timing of thermal conditions affected the vertical distribution and abundance of prey. Winter to early spring months in the cool year were characterized by a protracted winter stratification period (lasting until May as indicated by the onset of 4 °C isotherm), high March–April prey concentrations, and a later onset of summer stratification (Mid-May) (Figs. 10, 11). The warm year was characterized by an earlier break down of winter stratification and initiation of the spring isothermal period in April, an early onset of summer stratification in late April, warmer surface temperatures in Mid-May, and high prey concentrations in late April and May. Mussel grazing reduced the abundance of prey in the spring isothermal period when the whole water column at M45 was mixed. Summer stratification separated mussels from the SML thereby precluding them from filtering prey in the top 30 m (Figs. 10, 11).

Fig. 10
figure 10

Vertical distribution of temperature at the mid-depth location (M45) near Muskegon, MI from January-August for different meteorological years

Full size image
Fig. 11
figure 11

Vertical distribution of prey (phytoplankton, zooplankton, and detritus in J L−1) at the mid-depth location (M45) near Muskegon, MI for different meteorological years with and without mussels

Full size image

The cool year produced greater prey concentrations in March, early April, and early May than all other meteorological years, but cooler temperatures and mussel grazing rendered habitat unsuitable for both species in March and April (Figs. 12 and S2.9). At the mean April temperature in the cool year (2.7 ± 0.5 °C), Bighead and Silver Carp temperature-dependent filtration rates are 1564.2 L d−1 and 797.0 L d−1, respectively. At these diminished filtration rates, the model indicates that habitats must contain prey concentrations of 8.2 J L−1 to support Bighead Carp growth and 19.9 J L−1 for Silver Carp growth (the energetic equivalent of 3.15 μg L−1 and 7.7 μg L−1 Chla, respectively). Mean April prey biomass without mussels in the cool year (9.0 ± 1.4 J L−1) satisfied Bighead Carp’s growth requirement, but mussel filtration rendered the habitat unsuitable by reducing the mean prey biomass to 5.9 ± 2.2 J L−1. In May, when the model simulated the largest difference in Bighead Carp habitat suitability across the scenarios, the average temperature in the cool year in the top 30 m of the water column was 5.8 ± 0.96 °C. Temperature-dependent filtration rates at this temperature were 2481.7 L d−1 for Bighead Carp and 1313.7 L d−1 for Silver Carp, which translated to respective prey requirements of 5.9 J L−1 (2.3 μg L−1 Chla) and 14.9 J L−1 (5.7 μg L−1 Chla). Mean May prey biomass with mussels present in the cool scenario (6.8 ± 1.0 J L−1) was below Silver Carp’s energetic requirement, but not Bighead Carp’s—however, suitable Bighead Carp habitat was infrequent in May due to the lack of spatial and temporal overlap of sufficient prey biomass and water temperature.

Fig. 12
figure 12

Vertical distribution of Bighead Carp Hypophthalmichthys nobilis growth rate potential (GRP; g g−1d−1) at the mid-depth location (M45) near Muskegon, MI for different meteorological years with and without mussels

Full size image

In contrast, the warm year benefitted Bighead Carp GRP, but Silver Carp habitat was largely unsuitable throughout the year at M45 (Figs. 12 and S2.9). Average May temperature in the warm scenario in the top 30 m (8.0 ± 2.4 °C) effectively increased filtration capacity to 3391.1 L d−1 for Bighead Carp (36.6% increase from cool year FR) and 1845.4 L d−1 for Silver Carp (40.5% increase from cool year FR), which translated to a minimum prey requirement for growth of 4.9 J L−1 (1.9 μg L−1 Chla) and 12.5 J L−1 (4.8 μg L−1 Chla), respectively. Additionally, the earlier initiation of summer stratification in the warm year (May) mitigated the effect of mussel grazing by removing mussel access to the SML—thereby increasing mean prey biomass in the upper portion of the water column. The effect of SML separation from mussels on prey abundance increased from warm to cool scenarios. Mussel filtration accounted for a 3% reduction in average May prey biomass in the top 30 m of the water column in the warm scenario (without mussels = 7.5 ± 1.7 J L−1; with mussels = 7.3 ± 1.8 J L−1), a 10% reduction in the reference scenario (without mussels = 7.96 ± 1.3 J L−1; with mussels = 7.2 ± 1.3 J L−1), and a 25% reduction in the cool scenario (without mussels = 9.1 ± 1.4 J L−1; with mussels = 6.8 ± 1.0 J L−1). The timing of stratification simulated at M45 in the warm year coincided with the offshore expansion of suitable Bighead Carp habitat visualized in Fig. 6.

Discussion

Dreissenids: an impactful but imperfect control on BHC habitat

As we hypothesized, the effect of mussels on BHC habitat quality were temporally and spatially discrete—being limited to times and locations where they were in contact with the SML. The ability of mussels to clear the water column of phytoplankton and seston is a function of their biomass and depth, with the greatest clearance typically occurring in the 30–50 m depth range (Vanderploeg et al. 2010). Using biophysical models, Rowe et al. (2015, 2017) demonstrated how mussel filtration reduced phytoplankton biomass in the spring with the greatest impacts in the nearshore (15 m) and mid-depth areas (45 m), but had little effect during periods of stratification as impacts were reliant on the mussels being in contact with the SML. The consequent effect on BHC habitat quality that our model simulated reflects modeled and observed phytoplankton impacts (Vanderploeg et al. 2010; Rowe et al. 2015, 2017) and supports the notion that mussel grazing in the spring notably limits the quantity and quality of suitable BHC habitat, especially in the nearshore. However, the effect of food-limitation on BHC GRP due to mussel filtration is confounded by the effect of spring water temperatures. While cool temperatures decrease BHC respiration rates, which reduces the daily prey requirement needed for growth, they impose a relatively greater negative effect on foraging rate. Our model indicates that the reduction in foraging capacity (i.e. the fish feed less actively) due to colder temperatures requires BHC to occupy habitats with greater prey concentrations in a lake that is comparatively far less productive than what is typical of environments where they currently exist (Fahnenstiel et al. 2016; Tumolo and Flinn 2017; Deboer et al. 2018). This suggests that while mussels do reduce habitat quality in the spring, the degree to which they can effectively remove or limit suitable BHC habitat is contingent on other biophysical attributes (i.e. SML depth, water temperature, and total prey biomass) that are regulated by meteorology and nutrients.

Meteorological controls on thermal habitat and structure mitigate mussel impacts

The control of mixing dynamics on mussel impacts to BHC habitat suitability highlights how climate change can reduce the capacity of mussels to act as a barrier to BHC growth in certain areas of the lake. We found that a warmer winter-spring period enhanced habitat suitability by altering mixing dynamics and thermal characteristics to be more favorable for BHC growth. By advancing the onset of thermal stratification, a warmer winter-spring effectively reduces the duration of time BHC and dreissenid mussels would be in competition for food in the SML while simultaneously increasing BHC foraging rates due to increased temperatures. In our results, this effect was most apparent in deep water, where dreissenids’ access to the whole water column is restricted to periods of deep mixing (Rowe et al. 2015, 2017). The initiation of summer stratification at M45 that coincided with the timing of Bighead Carp habitat expansion into the offshore indicates that the earlier availability of suitable habitat in these regions was in part due to stratification isolating mussels from the SML. In fact, our results indicate that the effect of a warmer winter-spring on the average quantity and duration of suitable Bighead Carp habitat is comparable to what Lake Michigan’s suitability would be like if mussels were absent under the reference meteorological and nutrient conditions. The duration of summer stratification in Lake Michigan and other northern temperate lakes is expected to increase as climate warming progresses (Brooks and Zastrow 2002; Winder and Schindler 2004; Shimoda et al. 2011). Our findings indicate that this could further reduce the ability of mussels to act as a barrier to BHC establishment and effectively lengthen the period of time in which growth is possible.

While a warmer winter-spring improved habitat duration and quantity, Lake Michigan’s offshore waters would still limit Bighead Carp to near weight maintenance levels and would be incapable of supporting Silver Carp growth throughout the year. We found that a warmer climate led to greater accumulation of thermal energy in the winter-spring months, but mean annual CDD of Lake Michigan in a warm year (3813 ± 325 degree days in top 10 m) is still less than the current thermal regime in the western basin of Lake Erie (4433 degree days; (Kocovsky et al. 2012)), which is considered to be one of the most bioenergetically-suitable environments for BHC within the Great Lakes (Cooke and Hill 2010; Anderson et al. 2015). The warm year increased Bighead Carp growing days in less productive offshore waters, but both species experienced the greatest increase in annual growth in Green Bay—a comparatively shallower and more productive habitat than the rest of Lake Michigan. In Green Bay and other productive habitats, increased annual growth due to longer growing seasons could improve winter survivorship and fecundity (Degrandchamp et al. 2007; Coulter et al. 2018b), which is a dimension of establishment risk that we did not assess in this study. In the offshore, a warmer climate would not raise habitat quality into the same tier as the nearshore, but an increased duration of suitable habitat would create more opportunities for the spread and establishment of BHC. These fishes are capable of long distance migrations (Coulter et al. 2016) and can adjust to poor food conditions quickly—allowing them to fast for several weeks while minimizing weight loss (Sheng and Ma 2008; Coulter et al. 2018b). The longer availability of suitable offshore habitat, which would lessen the energetic burden of broad-scale movements, increases the likelihood BHC could survive migrations between productive habitats, which themselves would also be enhanced by a longer growing season.

The importance of nutrient loads to habitat suitability

While the invasion of the dreissenid mussels has altered productivity patterns and plankton abundance in Lake Michigan, nutrient load reductions have been the most impactful factor affecting BHC habitat quality over the past 50 years. The “High P” scenario is more typical of Lake Michigan’s total phosphorous loads in the 1970s and 1980s (Dolan and Chapra 2012) when efforts began to reduce total phosphorous loads in the Great Lakes. This scenario represents a 70% increase in loads from reference conditions, which as Rowe et al. (2017) determined, translated to a 24–36% increase in annual lake-wide mean primary production and increases of 0.8–1.1 μg L−1 in summer and fall surface chlorophyll. This increase in nutrients translated to Bighead Carp gaining an additional 2.2 kg and Silver Carp gaining an additional 0.64 kg over the course of year compared to the reference scenario. Even with mussels present, Bighead and Silver carp respectively gained 1.3 kg and 0.36 kg more weight over the year than in the reference scenario. Regardless of loading regime, loads had a greater effect than mussels on annual growth. This suggests that while the invasion of dreissenid mussels and nutrient load reductions have collectively reduced habitat quality over time, nutrient loads are the more influential factor and future nutrient enrichment would still improve suitable habitat quality despite the current mussel infestation.

Model-predicted annual growth in the nearshore of “High P” scenarios (Silver = 0.14–0.34 kg yr−1, Bighead = 3.0–3.8 kg yr−1 [range of mean nearshore annual GRP in the two High P scenarios]) is comparable to growth trends seen in habitats where BHC currently exist. Length-at-age data (converted to weights) from the Middle Mississippi River suggests that age 3 + Silver Carp are capable of growing 0.22–1.36 kg yr−1(Williamson and Garvey 2005) and age 3 + Bighead Carp tended to be 0.32–2.7 kg larger than the average weight of the previous year (Nuevo et al. 2004). While GRP models are intended to describe relative habitat quality instead of predicting realized growth (Tyler and Brandt 2001)—which requires consideration of other ecological factors such as habitat selection, inter- and intraspecific competition, and predation (e.g. Nelson et al. 2017; Coulter et al. 2018a, b)—contextualizing BHC habitat quality with reported growth rates can better communicate the implications of increased nutrient loads on Lake Michigan’s vulnerability to BHC. Furthermore, these results highlight the tight link between BHC habitat quality, primary production, and nutrient loads in phosphorus-limited systems like Lake Michigan (Shimoda et al. 2011; Warner and Lesht 2015; Rowe et al. 2017).

Tributary loads in Lake Michigan have consistently been under the 5600 MTA target and maintained at a steady level since the 1990s (Dolan and Chapra 2012), but changes in land use and river hydrology as a result of urban expansion and climate change could increase the amount of phosphorus loading into Lake Michigan. Wiley et al. (2010) found that average TP loads in the Muskegon River watershed could increase by 6.7–66.7% from 1998 levels by the end of the twenty-first century under scenarios of varying rates of urban expansion and climate warming. At a lake-wide scale, LaBeau et al. (2014) estimated that the TP load in Lake Michigan would increase by 4.9% from 2010 to 2040 based on expected urban expansion and by almost 10% when accounting for forecasted agricultural expansion in addition to urban sprawl. This increase in P loads is not substantial in the context of Lake Michigan’s 5600 MTA target, but it is significant to the establishment potential of BHC. Given that suitable habitat exists for both species under the current loading regime and remained extensive for Bighead Carp in our “No Loads” scenario, any future increases in TP loads would only create a more hospitable environment for BHC.

Model limitations

The model we used does not delineate phytoplankton or zooplankton by species composition or size class, nor did our foraging sub-model take into account size-selective BHC filtration and retention efficiencies. Differences in the size- and taxon-associated selectivity of dreissenids and BHC suggests that the degree of diet overlap between these invasive species would be more nuanced than our model assumes. Dreissenids are capable of filtering smaller particles than both species of BHC (Cremer and Smitherman 1980; Cotner et al. 1995; Tang et al. 2014), tend to be more efficient at reducing rotifers compared to larger zooplankton (MacIsaac et al. 1995; Jack and Thorp 2000; Higgins and Vander Zanden 2010), and are capable of promoting the dominance of toxic cyanobacteria such as Microcystis aeruginosa through selective rejection (Vanderploeg et al. 2001; De Stasio et al. 2014; Tang et al. 2014). In contrast, BHC have greater impacts on macrozooplankton resulting in the increase in relative abundance of rotifers (Sass et al. 2014), are capable of feeding throughout the water column regardless of the SML depth, and will consume Microcystis (Turker et al. 2003). As a result, both species may biomanipulate the phytoplankton and zooplankton communities in ways that further support their coexistence, which could lead to a collective ecological impact greater than mussels alone.

The biophysical model scenarios were developed to analyze the sensitivity of primary production to these stressors in Lake Michigan (Rowe et al. 2017), and therefore, are not representative of real lake conditions (except the reference scenario). Thus, our GRP model is not intended to be a true forecast, but rather a tool for understanding the implications of changes to Lake Michigan’s limnology for BHC establishment risk, and to inform and prioritize appropriate management actions. While the possible consequences of multiple anthropogenic stressors on BHC habitat quality are more numerous than our analysis could represent, our study provides a basis for inferring how changes in climate, land use, and the presence of dreissenid mussels could affect the risk of a BHC invasion in Lake Michigan.