© J.A. Dobosenski et al., Published by EDP Sciences, 2024

1 Introduction

Freshwater whitefish (Family Salmonidae; Subfamily Coregoninae) support economically and culturally important fisheries in northern North America and Eurasia (Bernatchez and Dodson, 1994; Anneville et al., 2015). This cold-water subfamily is particularly sensitive to overfishing, environmental degradation, climate change, and food web alterations stemming from invasive species, with perturbations resulting in a growing list of extinct, extirpated and threatened species (Vonlanthen et al., 2012; Anneville et al., 2015; Eshenroder et al., 2016). Like many species, survival of early coregonine life stages is paramount in setting cohort strength (Cushing, 1990; Eckmann, 2013). In some instances, management actions have mitigated anthropogenic stressors causing high mortality of early life stages. For example, eutrophication of some European lakes in the 20th century led to low dissolved oxygen concentrations on spawning grounds causing low survival of eggs and reductions in important egg incubation habitat (Müller, 1992; Eckmann, 2013) and ultimately loss of biodiversity (Vonlanthen et al., 2012). Modernization of sewage treatment in the latter half of the century led to re-oligotrophication of several lakes (Jenny et al., 2020) to the benefit of some coregonine populations (Anneville et al., 2009), but not all (Rellstab et al., 2004; Bourinet et al., 2023). Hatchery propagation programs were established to maintain fisheries in some European lakes with mixed success. Unfortunately, this intervention led to further loss of diversity in some lakes (Anneville et al., 2015). To stem the tide of global population decline, we need to better understand habitat requirements of early life stages of naturally reproducing populations (Eckmann, 2013) so that guidance for future conservation actions can reduce or eliminate unintended consequences.

The influence of abiotic and biotic factors on larval coregonine ecology has been studied for decades through a mix of field studies, mesocosm and aquarium experiments (Eckmann, 2013). Newly hatched larvae exhibit positive phototaxis and thermotaxis, with larvae of most species occupying surface waters during spring (Clady, 1976; Lahnsteiner and Wanzenböck, 2004; Myers et al., 2014). Moreover, larvae are visual predators (Mahjoub et al., 2008) with most feeding occurring during daytime and crepuscular periods (Johnson et al., 2009).

Lake trophic state can impact growth and survival of larvae as it influences production levels of early stages of zooplankton, especially copepod nauplii and copepodites that gape-limited hatching larvae use when exogenous feeding begins (Anneville et al., 2007, 2011; Davis and Todd, 1998; Lucke et al., 2020). In general, availability of prey is not thought to be a limiting factor to larval coregonines in most lakes but may be in systems which are ultra-oligotrophic (Eckmann, 2013) or significantly impacted by invasive filter-feeders (Cunningham and Dunlop, 2023).

Bathymetric depths used by larval coregonines can vary greatly across species and populations. In Lake Ontario (44.1490⁰, −77.2126⁰), cisco (Coregonus artedi) and lake whitefish (C. clupeaformis) larval densities were highest within 100 m of shore where bathymetric depths were <10 m (Brown et al., 2022). Similarly, European whitefish (C. lavaretus) larvae generally occur in the highest densities in lake littoral zones (Lahnsteiner and Wanzenböck, 2004). Conversely, cisco in Lake Superior (47.7231⁰, −86.9407⁰, Oyadomari and Auer, 2008; Myers et al., 2009; Lucke et al., 2020), alpine whitefish (C. zugensis) in Lake Lucerne (47.0136⁰, 8.4372⁰; Rellstab et al., 2004), and bloater (C. hoyi) in Lake Michigan (43.4501⁰, −87.222; Eppehimer et al., 2019) occupy greater distances from shore, and thus a comparatively wider range of bathymetric depths. Because larvae tend to be close to the surface (Myers et al., 2014), currents can transport larvae from spawning grounds (Hoagman, 1973), sometimes over great distances in large lakes (Oyadomari and Auer, 2008). Vendace (C. albula) in Lake Konnevesi, Finland, disperse eggs across the entire lake and the larvae continue to disperse after hatching with no correlation present between fall egg density and spring newly hatched larvae density among sites (Karjalainen et al., 2019, 2021). Another reoccurring theme is that larval densities can be highly variable, even at small spatial scales. This has befuddled researchers because the catch variation has not comported with obvious differences in habitat characteristics (Clady, 1976; Hoagman, 1973; Ryan and Crawford, 2014; Lucke et al., 2020). Water temperature has been associated with larval coregonid catches in some cases (Brown et al., 2022) but not others (Ventling-Schwank and Meng, 1995; Ryan and Crawford, 2014; Brown et al., 2022), despite the importance of temperature for growth, food availability, and mortality (Winder and Schindler, 2004; Eckmann, 2013; Ryan and Crawford, 2014; Jenny et al., 2020). Temperature was associated with larval lake whitefish catches but not with larval cisco in Lake Ontario (Brown et al., 2022), whereas no association was found between larval lake whitefish density and temperature in Lake Huron (Ryan and Crawford, 2014). Water temperatures can affect larval growth directly by regulating metabolic rates or indirectly by influencing behavioral performances (Eckmann, 2013). Lake temperatures influence the accrual rate of degree days of eggs which regulates hatch date and the size at hatching (Eckmann, 2013; Stewart et al., 2021). Water temperatures also influence food availability indirectly because during warm springs the phenology of plankton is accelerated (Winder and Schindler, 2004; Jenny et al., 2020). For example, earlier plankton blooms, stemming from an increasing trend in spring water temperatures, was tied to a substantial increase in European whitefish abundance in mesotrophic Lake Geneva (46.4414⁰, 6.5295⁰) during the 1990s (Anneville et al., 2009). Interannual variation in water temperatures had a strong influence on the sizes of European whitefish larvae at the end of spring in oligotrophic Lake Annecy (45.8165⁰, 6.1673⁰), with rapid cooling during spring stalling larval growth until warming resumed (Perrier et al., 2012). Rapid early spring warming can also lead to near complete mortality of lake whitefish larvae (Ryan and Crawford, 2014), whearas a cold winter can result in high overwinter survival of lake whitefish eggs, leading to intra-specific competition of larvae for food and decreased survival of this life stage (Freeberg et al., 1990).

Some version of the “bigger is better” paradigm is commonly invoked in the coregonine literature whereby conditions promoting larger hatching sizes and faster growth of larvae leads to earlier access to a broader range of prey sizes, increased swimming performance and burst speeds for prey capture and predator avoidance, and reduced predation risk as larvae more rapidly outgrow their predators’ gape limits (Rice et al., 1987; Eckmann, 2013). Consequently, one might expect the distribution of larvae should be strongly correlated with the density of edible-size zooplankton (EZ) prey. However, this has not always borne out in the literature. In support, the number of zooplankton in the stomachs of larval lake whitefish at a Lake Michigan site (43.2235⁰, 86.3365⁰) was correlated (R² = 0.51) with in situ measured cyclopoid densities (Pothoven, 2019). Conversely, larval cisco in Oneida Lake (43.2046⁰, 75.9231⁰) moved en masse from nearshore to offshore waters in spring despite the higher abundances of zooplankton prey in the vacated inshore waters (Clady, 1976). Moreover, zooplankton availability did not explain spatial variation in early summer catches of bloater larvae in Lake Michigan (Eppehimer et al., 2019).

Attempts have been made to relate larval growth, survival, and recruitment to the availability of zooplankton. For example, small changes in the number of copepods per larval lake whitefish in Lake Michigan resulted in large differences in instantaneous daily growth with negative growth occurring at <20 copepods per larvae and maximum asymptotic growth at >100 copepods per larvae (Freeberg et al., 1990). Growth of larval lake whitefish in the Bay of Quinte of Lake Ontario (44.1490⁰, 77.2126⁰) during the invasion of filter-feeding dreissenid mussels (Dreissena polymorpha and D. bugensis) was correlated with the number of EZ per larva (r = 0.813) and recruitment declined precipitously as edible cyclopoid copepod densities declined from a mean of 3.9 to 0.4 ind.L⁻¹ as invasive mussels expanded (Hoyle et al., 2011). Mortality of alpine whitefish larvae held in tanks for 34 days was low (<10%) when zooplankton densities were 150.ind.L⁻¹ but increased markedly to >35% when larvae were exposed to densities <20 ind.L⁻¹ (Rellstab et al., 2004). Re-oligotrophication of Lake Lucerne, Switzerland (47.014⁰, 8.437⁰) led to a decline in spring EZ concentrations from >20 ind.L⁻¹ in the 1970s to <10 ind.L⁻¹ after 1984, which coincided with a large drop in alpine whitefish year-class strength (Rellstab et al., 2004).

A multitude of factors likely influence larval distributions in complex ways, leading to many questions. Why can larvae be dispersed great distances through advection in some lakes (Oyadomari and Auer, 2004, 2008) whereas their distributions seem far from random with respect to bathymetric depth in others (Lahnsteiner and Wanzenböck, 2004; Brown et al., 2022, 2023)? Water temperatures have obvious bearing on larval food availability and foraging proficiency, but do larvae try and actively maintain association with favorable temperatures, or are they merely passive participants at the mercy of wind and currents? Availability of suitable-sized prey is important, especially in ultra-oligotrophic systems (Eckmann, 2013), but are larval densities correlated with available zooplankton densities at small spatial scales? To address these questions, we chose to study two lake systems that vary markedly in wind conditions, thermal regimes, and trophic state. We used repeated sampling at the time of hatching to determine whether and how larval distributions were related to habitat conditions. By studying two disparate systems in this way, we sought to (1) learn if abiotic and biotic variables shape larval distributions similarly and (2) use our results to attempt to answer these vexing questions.

2 Materials and methods

2.1 Study areas

Lake Geneva is the largest (580.1 km²) natural lake in western Europe (Fig. 1), with a mean depth of 153 m and an elevation of 374 m above sea level. Anthropogenic phosphorus inputs have been reduced since the 1970s, with in-lake total phosphorus concentrations declining from peak values of 80 μg L⁻¹ in the 1980s to 19 μg L⁻¹ by 2016 (Barbier et al., 2017) when we sampled.

Lake Superior is the world’s largest freshwater lake by surface area (82,100 km²), with an average depth of 147 m and an elevation of 180 m above sea level. Spring total phosphorous concentrations averaged 2.3 μg L⁻¹ between 1993 and 2010, making it ultra-oligotrophic (Barbiero et al., 2012).

For study sites, we chose a European whitefish spawning area near Thonon-les-Bains, France (Goulon et al., 2020) and a Lake Superior lake whitefish and cisco spawning area near Grand Portage, Minnesota, U.S.A. (Goodyear et al., 1982). Both study sites are characterized by limited littoral zones with a mix of boulder and cobble substrates and bathymetric depths >100 m within 2 km of shore.

Fig. 1 Locations of zooplankton (black circles) and larval tow stations (black lines) sampled during our Lakes Geneva (2016) and Superior (2018) campaigns. The Lake Geneva temperature logger at the laboratory port break wall is labeled INRAE (the National Research Institute for Agriculture, Food and the Environment laboratory), and the surface mini-logger station on Lake Superior is labeled Mini-logger.

2.2 Sample collections

At the Lake Geneva site, zooplankton and larval tows were conducted during six sampling events, spanning 4/4/2016 (month/day/year) to 4/21/2016; the Lake Superior site was sampled during eight weekly events from 5/3/2018 to 6/20/2018 once ice was out (Tab. 1). All sampling was conducted during the day, starting in the morning, and typically ending in the early afternoon (Tab. 1). We initially established two transects at Lake Geneva (one over 20–55 m bathymetric depths and a second over 60–100 m) and three transects at Lake Superior [15–50 m, 60–100 m, and 25–30 m (justification below)] to conduct larval fish tows (Fig. 1). Sampling initially included 8 and 10 larval tow stations at the Lake Geneva and Lake Superior sites, respectively (Fig. 1). The Lake Superior transects were placed outside Grand Portage Bay where a commercial fishery operates from June to November. Cisco egg densities in a nearby bay (Thunder Bay, Ontario; 48.495⁰, 89.06⁰) were highest over 20–30 m bathymetric depths having silt bottoms (Paufve et al., 2022) so we included the narrow 25–30 m transect range in the design. Zooplankton were collected at the midpoint of each larval tow (details follow).

In Lake Geneva, low larval catches led us to add a shallow-water zooplankton and larval tow station (Station 0, Fig. 1) on the 4th sampling event (4/14/2016) that we subsequently sampled. In Lake Superior, we used this same rationale to add three shallow larval tow stations (B1, B2, B3; Fig. 1) and a shallow zooplankton station (B1) on the 4th (5/24/2018) and subsequent sampling events. Sampled depths, after adding shallow sites, ranged from 6.8 m to 92.7 m on Lake Geneva, and from 3.9 to 93.6 m on Lake Superior (Tab. 2).

2.3 Surface temperatures

We collected either 4 or 8 water temperature profiles per sampling event on Lake Geneva, using either a multichannel logger (RBR Ltd., Ottawa, Ontario, Canada) or a Sonde EXO 1 probe (YSI Inc., Yellow Springs, Ohio, USA). Temperature profiles were collected either at all eight zooplankton stations (events 1 or 3) or two per transect (events 2, 4, 5 and 6) when time was limited. We supplemented these data with a profile at station 0 during event 6. We also obtained surface temperature data on an hourly timestep from a PT100 thermistor (Cimel Electronique, Paris, France) positioned on the National Research Institute for Agriculture, Food and the Environment (INRAE) laboratory port break wall (Fig. 1). For each probe deployment, we averaged all water temperature recordings in the uppermost 1 m of the water column where larvae were also collected.

The vessel used at our Lake Superior site lacked tackle to do overboard temperature profiles, so we deployed a HOBO® 8K pendant temperature data logger (Onset®, Bourne, Massachusetts, USA) at the surface at 47.942°, 89.658° (Fig. 1). This temperature logger recorded data every hour during the duration of our 2018 campaign. At the Lake Superior site, we calculated the mean surface temperature during each event based on the 24-hourly readings on the day of larval sampling.

2.4 Wind vectors

Wind direction and speed recordings were obtained from the nearest available weather station to each location. For Lake Geneva we used data obtained from Sciez weather station (46.331⁰, 6.3798⁰), roughly 10 km southwest of our study site (http://meteo-sciez.fr/site/index.php). For the Lake Superior site, we used the Rock of Ages (ROAM4) Coastal-Marine Automated Network (C-MAN) station, maintained by the National Data Buoy Center of the National Oceanic and Atmospheric Administration (NOAA). This weather station is located on the southwest end (47.867⁰, 89.313⁰) of Isle Royale, Michigan, USA, 29.5 km from Grand Portage, MN. Measurements of wind direction and wind speed were visualized by developing rose plots produced using the package ‘ggplot2’ (Wickham, 2009) in R (R Core Team, 2020). Data from two days prior to each sampling date and the first eight hours of the sampling day were used to build rose plots.

Cardinal wind directions in degrees offers some analytical challenges because the value defining a similar north wind direction can vary appreciably say from 10⁰ to 350⁰. We opted to use wind vectors (u for the east-west vector, and v for the north-south vector; Grange, 2014) that combine both wind direction and wind speed, using available measurements two days prior to each sampling date and the first eight hours of each sampling day. The choice of calculating wind vectors for 2.33 days prior to larval sampling events was admittedly arbitrary, so alternatives of averaging 1.33 and 0.33 days were also considered (see methods under Sensitivity of GAM modeling to explanatory data handling).

We calculated u and v vectors of each recording by first converting recorded wind direction from degrees to radians and then using the formulas: u = −Wind Speed (km.h⁻¹) * sin (Wind Direction), v = −Wind Speed (km.h⁻¹) * cos (Wind Direction). An average u and v were then calculated for each sampling event and these values were overlaid on the rose plots to aid interpretation. Zero wind speeds were excluded in the calculations of average u and v due to the zero wind speeds being assigned a default zero-degree direction in the datasets which would not be logical to include in the vector calculations. This equated to dropping 49% of wind observations for Lake Geneva only one wind observation for Lake Superior.

2.5 Larval sampling and processing

On Lake Geneva, we collected larvae with two ichthyoplankton nets, each having 1.5 × 1.0 m rectangular mouths, 5.0 m lengths, and 1-mm mesh (Colon et al., 2006; Anneville et al., 2007). To ensure the spacing and stability of the nets, two deflectors, at an angle of 30° with the front face of the frame, were positioned on the lower length and side of the frame. A collector was equipped at the end of the net and the mesh size was constant along the length of the net, ensuring the youngest larval stage of European Whitefish (size 10–12 mm total length, and >1 mm depth at hatching) were captured (Colon et al., 2006; Luczynski et al., 1988; Anneville et al., 2007). Nets were deployed with 4-m-long booms, so they fished the topmost 1 m of the water column on both sides of the vessel (midship) but outside the vessel’s wake.

We navigated between predetermined waypoints, targeting ∼15-minute tow durations and ∼400 m of tow distances. Effort was increased to 30-minutes of towing at station 0 during the 5th and 6th sampling event (∼800-m of effort) to increase catches to better characterize larval sizes. Vessel speed at Lake Geneva averaged 1.8 km.h⁻¹ (SD = 0.2 km.h⁻¹).

On Lake Superior, we collected larvae at Grand Portage with two conical nets, each having 0.5 m mouths, 1.5-m lengths, and 500-μm mesh (Myers et al., 2008; Brown et al., 2022). Nets were deployed with 1.0-m long booms so they fished the topmost 1 m of the water column. Deployment lines were 30-m long. Tows outside Grand Portage Bay (Fig. 1) were 5-minutes in duration and ∼350 meter of tow distance. The three tows added inside Grand Portage Bay were ∼10 min in duration and traversed ∼700 m of tow distance. Vessel speed at Grand Portage averaged 4.3 km.h⁻¹ (SD = 0.3 km.h⁻¹). Larvae from both lakes were placed in plastic bottles with lake water, labeled by station and vessel side, stored in coolers on ice, and transported to the lab for further processing, including transfer to 95% ethanol, within 24 h. Larvae were measured (when fresh) to the nearest mm and weighed to the nearest 0.001 g for Lake Geneva and 0.0001 g for Lake Superior. The Lake Geneva larvae were assigned to a larval developmental stage (LDS) using the algorithm developed for European whitefish: LDS = 0.3647 × length (mm) − 2.9342 (Luczynski et al., 1988).

Because different nets were used at the two lakes, we standardized larval densities by dividing larvae catches by the estimated volume of water sampled (ind.10³ m⁻³). Volume of water sampled was calculated by estimating distance traveled using GPS coordinates and assuming 100% net efficiency. Before pooling larval catches from both nets, we compared confidence intervals in mean catches between the port and starboard sides to determine if catch rates differed. Due to skewness of the data, confidence intervals were calculated using bootstrap resampling with 1000 replicates following the bias-corrected and accelerated (BCa) bootstrap method (Efron, 1987). The numbers of larvae caught in starboard- and port-side nets at the Lake Geneva site (starboard mean = 1.47, 95% CI = 0.94–2.35; port mean = 1.49, 95% CI = 0.94–2.69) and the Lake Superior site (starboard mean = 0.57, 95% CI = 0.13–2.56; port mean = 0.43; 95% CI = 0.16–1.31) did not differ, so samples were pooled for subsequent analyses. The high occurrence of zero catches, especially at the Lake Superior site, led us to add a constant, equal to one-half the smallest non-zero observed density at each lake (Leith et al., 2010). The smallest non-zero larval densities for Lake Geneva and Lake Superior were 0.6 ind.10³ m⁻³ and 3.2 × ind.10³ m⁻³, respectively, so we added 0.3 to zero larval catches on Lake Geneva and 1.6 to Lake Superior prior to data visualization.

We calculated mean and median larval densities with 95% confidence intervals and lower (25%) and upper quartile (75%) ranges across all samples and for three bathymetric depth ranges [<10 m (both lakes), 20–55 m (Lake Geneva) and 16–55 m (Lake Superior), and 60–100 m (both lakes)] by sampling event and for all sampling events pooled. This bathymetric range binning allowed us to explore how larval densities varied by bathymetric depth. Due to the skewness of the data, we calculated confidence intervals of the mean by using the BCa method (Efron, 1987) with 1000 replicates for all except the shallow bin due to small sample size.

2.6 Species identification of Lake Superior larvae

Mouth gape at the larval stage varies across predominant coregonine species found in Lake Superior (Davis and Todd, 1998) so identifying larvae to species was important. Larvae collected inside Grand Portage Bay during the 7th sampling event (6/14/2018, n = 37) were assigned to species using a Genotyping‐in‐Thousands by sequencing (GTseq; Campbell et al., 2015) panel at the U.S. Geological Survey’s Great Lakes Science Center, Molecular Ecology Laboratory (Ann Arbor, Michigan, USA) following the protocols outlined in Weidel et al. (2022).

2.7 Zooplankton sampling and processing

Zooplankton samples at Lake Geneva during event 1 (4/4/2016) were collected with a 64-μm closeable zooplankton net having a 0.5-m diameter mouth and a length of 2.5 m. The net was unmanageable given its length, so during event 2, we collected zooplankton with the closeable net and the standard alpine-lake Bongo net used in long-term surveys (Rimet et al., 2020) to develop a correction factor. The Bongo double net had mesh sizes of 64 μm and 200 μm, 0.35 m mouths and a length of 1.0 m. During all subsequent events of Lake Geneva and all eight sampling events of Lake Superior, we employed the alpine-lake Bongo net. Standard tows from 25 m to surface, or from 1.5 to 2.5 m above the lakebed to the surface at shallower sites, were conducted at each zooplankton station during each sampling event with the following exceptions. The net comparison during event 2 at Lake Geneva limited collections to stations 2, 4, 5 and 8, and high wind prevented sampling of sites D7, D8 and D9 during the 3th Lake Superior sampling event. At both lakes, zooplankton nets were washed into the net cod ends and samples were preserved in 95% ethanol, labeled, and stored in plastic bottles.

To understand the proportion of zooplankton residing in the uppermost water column where larvae were collected, we collected some shallow tows from 4 m to the surface (N = 2) and 6 m to the surface (N = 1) on Lake Geneva. A total of 14 shallow tows from 5 m to the surface were collected on Lake Superior. Shallow tows were collected right after the standard 25 m to surface tows.

Zooplankton samples were counted and measured at the Rubenstein Ecosystem Science Laboratory (University of Vermont, Burlington, Vermont, USA). Water was added to each sample to a known volume, from which 5–18 ml aliquots were drawn randomly for identification, counts and length measurements. Counts from the aliquots were used to estimate numbers in the sample, and then scaled to density (ind.L⁻¹) based on net sample volume assuming a filtering efficiency of 1. Zooplankton were measured under an Olympus SZX12 dissecting scope (Shinjuku City, Tokyo, Japan) equipped with a drawing attachment and digitizing tablet. Adult copepods and cladocerans were classified to genus, immature copepods were classified as calanoid or cyclopoid copepodites, and nauplii were assigned as their own group. Nauplii were not measured in the Lake Superior samples.

During quality assurance procedures, we identified that one sample processor had measured lengths to the end of the caudal setae (copepods) and tail spines (cladocerans) for a subset of samples from both lakes. To correct this error, we developed correction factors (see Appendix A for full details) to convert total lengths (to end of setae or tail spine) to standard lengths (to end of caudal rami or base of the tail spine).

2.8 Definition of edible zooplankton (EZ)

Because larvae are gape-limited, we needed to define for each lake the zooplankton species and sizes of prey that were edible. We considered the larval species, expected larvae consume prey slightly smaller prey than their mouth gapes (Schael et al., 1991), conducted a literature review of taxa of zooplankton consumed by larvae in comparable alpine systems to Lake Geneva (Anneville et al., 2007, 2011), and used results from a recent Lake Superior larvae diet study (Lucke et al., 2020). The methods to define EZ are provided in Appendix B. We opted to use results of the 64-μm-mesh net tows when estimating EZ densities because of the smaller size of zooplankton this mesh captures compared to 200-μm mesh. In Lake Geneva, EZ were defined as cyclopoid copepods ≤0.7 mm, cladocerans ≤0.6 mm, and nauplii. Calanoid copepodites and adults, regardless of size, were excluded (see Appendix B for details). In Lake Superior, EZ were defined as copepod nauplii, calanoid copepodites and adults < 0.75 mm. Cladocerans and cyclopoid copepodites and adults, regardless of size, were excluded (Appendix B).

Once defined, EZ densities were calculated for each sample. We also calculated mean and median EZ densities with confidence intervals and lower (25%) and upper quartile (75%) ranges for three bathymetric depth ranges [<10 m (both lakes), 20–55 m (Lake Geneva) and 16–55 m (Lake Superior), and 60–100 m (both lakes)] by sampling event and for all sampling events pooled. This bathymetric range binning allowed us to explore how EZ varied by bathymetric depth using larger sample sizes. Due to the skewness of the data, we calculated confidence intervals of the mean by using the BCa method (Efron, 1987) with 1000 replicates for all except the shallow bins due to small sample size.

2.9 Analysis of abiotic and biotic factors influencing larval spatial distributions

Generalized additive models (GAMs) were used to explore abiotic and biotic factors shaping larval coregonine distributions in both lakes. GAMs allow for linear and nonlinear relationships between individual explanatory variables and a response variable using smoothing functions. Due to the non-linearity within our data, we chose to use GAMs over generalized linear models (GLM). We used the R package ‘mgcv’ to develop and assess various GAMs for both lakes (Wood, 2015). The R package ‘fitdistrplus’ was used to examine the distribution of the response variable, larval coregonine count, against fitted distribution functions (e.g. Poisson, negative binomial, etc.) to determine the most appropriate family distribution (Delignette-Muller and Dutang, 2015). Volume was used as an offset in each model which enabled us to directly model count data but adjust for variation in total volume sampled per sampling event. Gam.check was used to produce diagnostic plots that identified potential issues with model fitting when testing multiple family distributions and link functions for each lake. Gam.check was also used to optimize smoothness selection (k) for each predictor (Wood, 2006). Correlation plots revealed that the u and v wind vectors exhibited collinearity (data not shown), necessitating they be loaded in separate models. Models were developed to encompass all possible combinations of predictor variables and then model.sel of the R package ‘MuMin’ was used to compare all models (Anderson and Burnham, 2002; Burnham and Anderson, 2004). Model.sel identified the most parsimonious model(s) by ranking all models based on Akaike information criterion corrected for small sample size (AICc). This method for identifying the best models followed that of a similar study by Eppehimer et al. (2019). We used the closest available EZ estimate when a given larval catch lacked paired measurements. For example, estimates of EZ measured at B1 at Grand Portage were applied to the B2 and B3 larval catches where no zooplankton had been collected (Tab. 2). Our full models for predicting larval counts in both lakes included EZ densities, bathymetric depth, Julian Day, u and v. A smoothing term function with a smooth class of thin plate regression splines was applied for all predictor variables. Surface temperature was also considered as a possible predictor but due to the high correlation between temperature and Julian day, we could not include surface temperatures in the models. However, we investigated the effect of temperature on larval density separately for Lake Geneva by using the closest available information (water column probe or INRAE thermistor) to each larval tow on each date to assess the role of temperature on larval densities in Lake Geneva.

For Lake Geneva the full models were:

Larval Coregonid Count ∼ s(EZ Density, k=10)+s(Bathymetric Depth, k=9)+s(Julian Day, k=6)+s(v, k=6), family=negative binomial (link=“log”), offset=log(Volume);

Larval Coregonid Count ∼ s(EZ Density, k=10)+s(Bathymetric Depth, k=9)+s(Julian Day, k=6)+s(u, k=6), family= negative binomial (link=“log”), offset=log(Volume).

For Lake Superior the full models were:

Larval Coregonid Count ∼ s(EZ Density, k=10)+s(Bathymetric Depth, k=10)+s(Julian Day, k=8)+s(v, k=8), family= negative binomial (link=“log”), offset=log(Volume);

Larval Coregonid Count ∼ s(EZ Density, k=10)+s(Bathymetric Depth, k=10)+s(Julian Day, k=6)+s(u, k=8), family= negative binomial (link=“log”), offset=log(Volume).

We report the approximate significance of smooth terms, indexed by P-values derived from Wald tests. Only models with ΔAICc< 4 are provided for Lake Superior and ΔAICc < 1 for Lake Geneva. We visualized the top-performing model in each lake with biplots of larval density against predictor variables included in said model.

2.10 Sensitivity of GAMs to explanatory data handling

When performing the GAM modeling to explain variation in larval distributions, we made several decisions about the handling of explanatory variables that could influence outcomes. First, we occasionally borrowed EZ densities from nearest neighbors when data were unavailable for Lake Superior, so we repeated the GAMs using only samples with no borrowing. Second, we explored how sensitive results were to our choice of averaging wind vectors for 2.33 days prior to larval sampling by building GAMs with wind vectors averaged for 1.33 days prior to sampling on both lakes, and for 0.33 days on Lake Superior only. Because station wind measurements in the 8 hours preceding the 1st Lake Geneva sampling event were all zero, the wind vectors could not be calculated. Additionally, we wanted to investigate how dropping the zero wind values might impact our results for Lake Geneva given the high percent of zero values (49%). This was not needed for Lake Superior given only one zero observation was present in the dataset. To investigate this, we created models with mean wind speeds (including zero values) per sampling event and investigated the importance of wind speed when all values are included. We also recognized that, the 64-μm mesh zooplankton net could be less efficient than the 200-μm mesh so we built GAMs using the 200 μm EZ density estimates while keeping our definitions of EZ the same. A lack of 200-μm net sampling during the 1st sampling event on Lake Geneva necessitated dropping this event from this alternative. Finally, we acknowledge that the transects for the three depths zones at each lake were close together. Therefore, we decided to run the models with values averaged for each depth zone since the transects proximity to each other might violate true independence.

3 Results

3.1 Surface temperatures

Lake Geneva had considerably warmer spring surface water temperatures because it is at a lower latitude, does not freeze over during winter, and has a smaller surface area than Lake Superior (Fig. 2). Surface water temperatures at Lake Geneva during April 2016 ranged from 8.9 to 12.8 °C. Surface temperatures at the Lake Superior site increased gradually from an average of 1.7 °C on 5/3/2018 to 3.6 °C on 6/20/2018 (Fig. 2). Surface temperature and date were correlated for Lake Geneva (r(30) = 0.82, p < 0.001) and Lake Superior (r(6) = 0.99, p < 0.001).

Fig. 2 Surface temperatures (Celsius, °C) versus calendar date (MM/DD) at Lake Geneva sites sampled during 2016 and Lake Superior sites during 2018. The Lake Geneva temperatures are the averages from 0–1 m depths from water column profiles and the thermistor positioned at 1m depth at the National Research Institute for Agriculture, Food and the Environment laboratory (INRAE) laboratory break wall (Fig. 1). The Lake Superior temperatures are the average of 24-hourly recordings measured by a data logger deployed at the surface at a single site (Mini-logger; Fig. 1).

Fig. 3 Larval coregonid natural log-transformed densities (ind.103m−3) and total lengths (mm) for each larval fish by sampling date in Lake Geneva and Lake Superior. Shallow site (<10 m bathymetric depths) larval densities (Lake Geneva site 0 and Lake Superior sites B1, B2 and B3, Fig. 1) and lengths are shown with open diamonds. Before transforming the data, half of the smallest non-zero value were added for each lake (0.3 for Lake Geneva and 1.6 for Lake Superior) due to the presence of zeros (Leith et al., 2010) which places zero values at −1.20 for Lake Geneva and 0.47 for Lake Superior after transformation. The black line within the box marks the median, the boundary of the box closest to zero indicates the 25th percentile, and the boundary of the box farthest from zero indicates the 75th percentile. Whiskers represent 1.5 * IQR from the bottom and top of the boxes, where IQR is the inter-quartile range or distance between the first and third quartiles. Data beyond the whiskers are outliers.

3.2 Wind vectors

Wind speeds for Lake Geneva ranged from 0 to 20.9 km.h⁻¹ with variable wind directions (Supplemental Fig. 1A). Minimal wind occurred in the 2.33 days leading up to the first three sampling events of Lake Geneva. Wind direction was dominant from the south prior to the 4/14/2016 and 4/18/2016 sampling events and was dominant from the north prior to the 4/21/16 sampling event. The highest winds were measured prior to the 4/18/2016 sampling event and blew from the coast towards the lake.

Wind speeds for Lake Superior ranged from 0 to 71.3 km.h⁻¹ with variable wind directions (Supplemental Fig. 1B). Wind direction was dominant from the northeast prior to the 5/9/2018, 5/18/2018 and 5/31/2018 sampling events, southwest prior to the 5/24/2018 and 6/14/2018 sampling events, and variable and light prior to the 5/3/2018, 6/7/2018 and 6/20/2018 sampling events.

3.3 Larval catches and sizes

Lake Superior had a much higher percentage of zero-catch tows (82.1%) than Lake Geneva (25.5%). A total of 148 larvae were collected at the Lake Geneva site in 51 tows (after pooling the two sides). Average density was 1.96 ind.10³ m⁻³ (95% CI = 1.38–2.80 ind.10³ m⁻³) over all sampling events with a median value of 0.87 ind.10³ m⁻³ (Tab. 3). Average larval densities generally trended upward from the 4/4/2016 to 4/18/2016 sampling events but declined sharply on 4/21/2016 (Fig. 3). Total length was measured for 143 of the larvae that were sufficiently intact. The sizes of larvae caught at the shallowest site sampled (site 0) after 4/14/2016 were comparable to sizes caught at the deeper sites during this timeframe (Fig. 3). Most of the measured larvae (94%) were 10-20 mm during all events (Fig. 3) with lengths ranging from 11.5 to 25.0 mm, averaging 15.9 mm over all sampling events. 81% of larval fish had predicted LDS scores of 1–3. The threshold for LDS values less than or equal to 3 is demarcated by the start of the development of the fin rays (Luczynski et al., 1988).

A total of 95 larvae were collected at the Lake Superior site in 95 tows; 86 larvae (91%) were caught inside Grand Portage Bay (sites B1, B2 and B3, Figs. 1 and 3) where bathymetric depths were < 10 m. Average larval density across all samples was 3.94 ind.10³ m⁻³ (95% CI = 1.45–13.27 ind.10³ m⁻³) with a median value of 0.00 ind.10³ m⁻³ due to the high numbers of zeroes (Tab. 3). Average larval density over all sampling events in the bay was 20.29 ind.10³ m⁻³ (median = 3.53 ind.10³ m⁻³), whearas densities for outside the bay (stations S1-S2, M3-M6, D7-D10, Fig. 1) averaged 0.68 ind.10³ m⁻³ (median = 0.00 ind.10³ m⁻³). Most larvae caught at the Lake Superior site were ≤14 mm (88%), with all larvae captured ranging from 9.0 to 22.0 mm (Fig. 3) and averaging 12.9 mm.

3.4 Lake Superior larvae species assignments

Of the 37 larvae caught on 6/14/2018 in Grand Portage Bay, 36 were assigned as cisco (97.3%) and one as kiyi (C. kiyi, 2.7%) via GT-seq with all assignment probabilities ≥0.965. The mean length of the cisco larvae was 13.2 mm with a range of 11–15 mm and the one kiyi larva had a length of 12 mm. This finding is consistent with knowledge that Grand Portage area is a cisco spawning area (Goodyear et al., 1982).

3.5 Edible zooplankton densities

A total of 52 and 79 64-μm mesh zooplankton samples were processed from Lakes Geneva and Superior, respectively (Tab. 2). Linear models developed to predict standard lengths of zooplankton from total lengths all had R² >0.77 (Appendix A, Tab. A1), and these models were applied to the samples with incorrectly measured total length. The conversion factor used to correct densities collected with the large closeable net on 4/4/2016 on Lake Geneva was density (ind .L⁻¹) = 0.01988 + 0.46886 (Closeable Net Density, ind.L⁻¹), R² = 0.92, F(1,14) = 164.1, p < 0.0001).

The average EZ density over all Lake Geneva sites and events was 4.42 ind.L⁻¹. The EZ was dominated by nauplii, with generally lower densities of cyclopoid copepods and cladocerans (Fig. 4). Average density of EZ for Lake Geneva larvae was highest at <10 m bathymetric depths (mean = 7.12 ind.L⁻¹, Tab. 3) and ranged from 2.78 to 11.90 ind.L⁻¹ during the last three sampling events (4/14/2016–4/21/2016) when densities were measured at shallow depths. During the entire April campaign, EZ densities for larvae were similar at the 20–55 m (mean = 3.87 ind.L⁻¹, range = 1.80 to 6.94 ind.L⁻¹ (Tab. 3)) and 60-100 m bathymetric strata (mean = 4.66 ind.L⁻¹, range = 1.17 to 8.21 ind.L⁻¹ (Tab. 3)).

The EZ densities for cisco larvae were considerably lower at the Lake Superior study area (Fig. 5) compared to Lake Geneva and averaged 1.4 × 10⁻² ind.L⁻¹ over the 2018 campaign (Tab. 3). The relative densities of copepod nauplii and small calanoid copepods varied during the Lake Superior campaign with no systematic pattern (Fig. 5). The EZ densities did not vary appreciably across the three bathymetric ranges.

Analysis of 14 paired Lake Superior tows from 25 m and 5 m to the surface, respectively, showed EZ densities were similar across these two depths (data not shown). A paired t-test showed the average density of EZ in our standard 25-m vertical tows (1.2 × 10⁻² ind.L⁻¹) and 5-m vertical tows (1.3 × 10⁻² ind.L⁻¹) did not vary significantly (t = –0.322, df = 13, P = 0.75), indicating densities of EZ did not vary between these two depth layers. We deemed the sample size of surface tows in Lake Geneva (N = 3) too small to warrant statistical comparison.

Fig. 4 Edible zooplankton densities (ind .L−1) for larvae based on 64-μm mesh tows at the Lake Geneva stations as a function of bathymetric depth range and sampling date. Stations were sampled between 4/4/2016 and 4/21/2016. Stations were binned as follows: <10 m bathymetric depths = site 0; 20–55 m bathymetric depths = sites 1–4; 60–100 m bathymetric depths = sites 5–8 (see Fig. 1). Stations <10 m were added on 4/14/2016. The 4/4/16 samples were collected using a closeable net instead of a double Bongo net, so density estimates were adjusted based on a conversion factor we developed (see Results section). Edible zooplankton on Lake Geneva was defined as cyclopoid copepods ≤0.7 mm, cladocerans ≤0.6 mm, and nauplii (see Appendix B for details).

Fig. 5 Edible zooplankton densities (ind L−1) for larvae based on 64-μm mesh tows at the Lake Superior stations by bathymetric depth range and sampling date. Stations were sampled between 5/3/2018 and 6/20/2018. Stations were binned as follows: < 10 m bathymetric depths = B1; 20-55 m bathymetric depths = sites S1-S2 and M3-M6; 60-100 m bathymetric depths = sites D7-D10 (see Fig. 1). The < 10 m zooplankton station was added on 5/31/2018. Edible zooplankton on Lake Superior was defined as copepod nauplii, calanoid copepodites and adults <0.75 mm; cladocerans and cyclopoid copepods, regardless of size, were excluded (see Appendix B for details).

Fig. 6 Generalized additive model predicting surface larval counts in Lake Geneva (left panel) and larval counts with volume as an offset in Lake Superior (right panel) as a function of explanatory variables. Circles are larval density observations and shaded areas are 95% confidence intervals. Jittering was applied when plotting due to overlapping data points. A break in the y-axis was applied for Lake Superior due to a high larval density observation and horizontal lines are provided to accentuate the break.

3.6 Analysis of abiotic and biotic factors influencing larval spatial distributions

For both Lake Superior and Lake Geneva, a negative binomial family was the best fit for the larval coregonid count data (Supplemental Fig. 2). For Lake Geneva, the negative binomial family with a log link led to the best diagnostic plots with the quantile-quantile (QQ) plot and histogram showing general normality of residuals, no pattern between the residuals and the linear predictor, and close to a 1:1 pattern between the response and fitted values. For Lake Superior, a negative binomial family with a log link led to a normal distribution of residuals based on the QQ plot and histogram, but the residuals vs. linear predictor plot indicated a pattern and a 1:1 relationship for all values in the response vs. fitted values plot was not evident.

Five Lake Geneva GAMs had ΔAICc within 1 of the top-performing models (Tab. 4). Of the six top-performing Lake Geneva models, all included bathymetric depth, whereas Julian Day, u (east-west wind vector) and v (north-south wind vector) were present in three of the models, and EZ density was present in two models. The weight of the top-performing model (0.164) was close to the weight of the 2nd model (0.144). Approximate P-values of the smooth terms of the top-performing Lake Geneva models were significant for bathymetric depth (<0.001) and v (<0.01) but not for Julian Day (0.09), which provides strong evidence for inclusion of bathymetric depth and u but less support for the inclusion of Julian Day. Based on the top-performing model, larval counts in Lake Geneva declined with bathymetric depth in a negative exponential fashion, and with a slight linear decline as Julian Day increased (Fig. 6). There was also an increase in larval count as v becomes positive which is related to southerly winds (Fig. 6). EZ density and u were not in our top model for Lake Geneva (Fig. 7).

The four top-performing Lake Superior models included bathymetric depth and u four times, and Julian Day and EZ twice (Tab. 4). The weight of the top-performing model (0.511) was higher than the 2nd model (0.200). Based on the top-performing model, larval counts dropped sharply with bathymetric depth with predicted counts near zero beyond ∼10 m depth (Fig. 6). To ensure that the highest larval count was not solely driving the GAMs results, we ran the model with the highest larval count (59) removed. After dropping this value, the top model stayed the same and depth and u were still significant in the top performing models. Some of the highest catches of larvae at the Lake Superior site occurred when u was >10 (Fig. 6), coinciding with a light southwest wind prior to the 6/14/2018 sampling event (Supplemental Fig. 1B) with the wind blowing towards the coast. Larval counts were also elevated at sites outside Grand Portage Bay (Fig. 3) when u was <9 (Fig. 6) during the 5/9/2016, 5/18/2016 and 5/31/2016 sampling events when winds were predominantly from the northwest blowing away from the coast. Larval counts at our Lake Superior site were unrelated to Julian Day, v or EZ densities (Fig. 7).

Fig. 7 Observations for predictor variables not included in the top models for each lake vs. larval density. Jittering was applied when plotting due to overlapping data points.

3.6.1 Sensitivity of GAMs to explanatory data handling

In general, the Lake Geneva GAM modeling was insensitive to the wind averaging decision (Tab. 6). When Lake Geneva wind vectors were averaged for 1.33 days preceding larval collections, the same three explanatory variables (bathymetric depth, Julian Day, and v) were identified as important in this alternative model (Tab. 5). For Lake Geneva, when we modeled wind speed instead of wind vectors to investigate the impact of dropping zero wind speed values, depth and wind speed were important while Julian Day was less important with higher wind speeds and lower depths being related to the highest larval catches (Tab. 5). When EZ densities were calculated using the 200-μm mesh, the GAMs showed strong support for including u (east-west wind vector), whereas Julian Day and v became less important (Tab. 5). When average values were used for shallow, mid, and deep transects, depth and v were identified as the most important explanatory variables (Tab. 5).

Beyond the GAMs described in our methods, we ran 4 additional alternative Lake Superior models (Tab. 5). When we did not borrow missing explanatory variables, the GAMs produced the same results. Using the 200-μm mesh for EZ or averaging wind vectors for 1.33 days and 0.33 days prior to larval sampling also had no impact on the top performing model (Tab. 5). When average values were used for shallow, mid, and deep transects, depth was identified as the most important explanatory variable (Tab. 5).

4 Discussion

Our analysis of factors influencing larval coregonine distributions revealed that bathymetric depth was important in both systems − the highest densities were observed where depths were <10 m (Fig. 6). Julian Day helped explain larval distributions to some extent in Lake Geneva, but not in Lake Superior. The v (north-south) wind vector was an important driver of larval distribution in Lake Geneva but the u (east-west) vector was important at our Lake Superior site. Densities of EZ did not seem to be an important driver of larval distributions in either system. Beyond the importance of bathymetric depth in both systems, the importance of the other abiotic and biotic factors we explored were inconsistent.

Our findings are consistent with prior European whitefish studies (Anneville et al., 2007, 2011; Lahnsteiner and Wanzenböck, 2004) and Lake Ontario cisco studies that identified shallow water as critical larval nursery habitat (McKenna et al., 2020; Brown et al., 2022, 2023). The high larval cisco densities we observed at our shallowest Grand Portage site stands in sharp contrast with prior Lake Superior studies that showed larval cisco were more abundant over greater bathymetric depths (Oyadomari and Auer, 2008; Myers et al., 2009, 2014; Lucke et al., 2020), however, only Lucke et al. (2020) reported sampling at depths <10 m and Oyadomari and Auer (2008) did not report depths. A recent analysis of genetic structure in Laurentian Great Lakes cisco using microsatellites concluded that cisco in western Lake Superior were likely one large metapopulation (Stott et al., 2022). Of the nine Lake Superior sites that were included, however, Grand Portage was identified as one of the most unique due to significant estimates of pair-wise differentiation or fixation index (FST) with three sites, compared to the apparent panmixia of the other eight sites evaluated (only one other Lake Superior site had a single, significant pair-wise FST estimate). Stott et al. (2022) hypothesized that a general absence of stock structure in Lake Superior cisco compared to Lake Huron, where stock structure was evident, could be tied to differences in the ability of larvae in these two systems to maintain position near hatching grounds. This is consistent with a recent study showing unfertilized cisco eggs were significantly larger in Lakes Huron and Ontario compared to Lake Superior (Koenigbauer et al., 2022), assuming larger eggs equate to larger more powerful larvae. A similar hypothesis has also been formulated for Atlantic herring (Clupea harengus) where the number of distinct stocks is determined by larval retention zones that are possibly delineated by hydrographic features (Iles and Sinclair, 1982). Whether the significant genetic differences between Grand Portage and several other Lake Superior sites detected by Stott et al. (2022) could be related to the unusually shallow bathymetric distribution of Grand Portage cisco larvae we collected in 2018 would need to be further assessed over broader spatial and temporal scales. Clearly, shallow water represents critical nursery habitat for larval coregonines in Lake Geneva, and in the case of Grand Portage, this habitat may be more important than past Lake Superior cisco studies have indicated.

Julian Day explained some of the variation in larval distributions at our Lake Geneva site, but not at our Lake Superior site. Based on the model, larval counts in Lake Geneva surface waters trended downwards throughout the campaign (Fig. 6). One explanation of this could be that larvae are dying throughout the time period of our sampling causing the decline in density. This was the first study investigating the distribution of coregonine larvae in Lake Geneva. Therefore, we do not have knowledge of typical temporal within-season dynamics of larval density for Lake Geneva to compare our observation of a decline in larval density throughout the spring to. In nearby Lake Annecy, there was no clear pattern in larval abundance throughout spring (Anneville et al., 2011) however, the two systems are not very similar in terms of surface, depth, and wind conditions. For Lake Superior, we do not have much knowledge of species specific within-season dynamics of larval densities due to the difficulties in identification larval coregonines to the species level at the larval stage. Recent advances in genetic identifications of these species have enabled us to begin to evaluate species specific information such as hatch time, habitat occupation, growth, and mortality rates (Lachance et al., 2021). In the Apostle Islands region of Lake Superior, larval cisco appeared during the week of May 21 and declined after the week of June 18 with length-frequency distributions indicating that they generally grew as a cohort until early July (Lachance et al., 2021). Julian day was not identified as an explanatory variable of larval density for our study, however, we didn’t sample as late into the summer as Lachance et al. (2021).

Another explanation could be related to temperature. Given that temperature data was strongly correlated with Julian Day (r(30) = 0.82, p < 0.001), we considered that temperature might be an alternative explanation for this trend. The highest catches occurred when surface temperatures were between 9° and 12 °C, a pattern consistent with previous observations from Lake Annecy, France (Perrier et al., 2012). Metabolic rates increase with temperature and require increased prey consumption to achieve growth (Karjalainen, 1992; Myers et al., 2014). Consequently, we postulate larvae in Lake Geneva may have switched to deeper habitat and cooler temperatures at the end of our 2016 campaign. We were unable to test if surface temperatures in Grand Portage Bay, where larval catches were highest, were different than outside the bay.

Because food availability is expected to be more critical to larval survival in ultra-oligotrophic lakes (Eckmann, 2013), we expected that larvae at our Lake Superior site may attempt to reside in areas with the highest EZ densities to avoid starvation. Such behavior may be less critical in Lake Geneva given their prey is generally not in short supply (Fig. 4). When we analyzed data with the GAMs, each combination of date and site was considered an observation and larval counts were unrelated to EZ at Grand Portage (Fig. 7). The EZ density was included in our fifth and sixth Lake Geneva GAMs (Tab. 4) with it being significant in one of them. The highest larvae counts occurred where EZ densities ranged from 2.5 to 7.5 ind.L⁻¹ (Fig. 6). However, in both systems, bin averaging by bathymetric stratum showed EZ densities were highest in the shallowest stratum (<10 m bathymetric depths) studied, and, interestingly, these shallow depths also supported the highest larval densities on average (Tab. 3). Admittedly, the mean EZ in the <10 stratum at Grand Portage was based on a low sample size (N = 3) and was driven by an exceptionally high value obtained on the last event. One possible explanation for our findings is that larvae do try to maintain position in shallow bathymetric depths to generally increase their feeding success but are not able to associate with patches of high EZ available within this shallow water foraging arena. Given we measured EZ at fixed sites and collected larvae over long tow distances, we cannot say definitively that larvae were not associated with patches of high EZ availability, but we consider this highly unlikely given our intense zooplankton sampling did not reveal the existence of such patches.

Wind was an important driver of larval distributions in Lake Geneva and Lake Superior. Larval densities at our Lake Geneva site peaked in the shallowest depths during our 4/18/2016 sampling event. Interestingly, this event was proceeded by a strong southerly wind which was the highest measured during our campaign that would have acted to blow surface waters offshore (Supplementary Fig. 1A). Despite this strong wind event, larval European whitefish were still found in the highest densities at shallow bathymetric depths on 4/18/2016. Conversely, cisco larval densities at our Lake Superior site were highest at the shallowest depths we sampled inside Grand Portage Bay during the 6/14/2018 sampling event when moderate southwest winds blew towards the coast prior to sampling. Prior to that, some strong winds blowing offshore from the northwest during the 5/9/18, 5/18/2018 and 5/31/2018 sampling events (Supplemental Fig. 1B) coincided with slightly elevated densities at the deepest most offshore sites we sampled (Figs. 3 and 6).

So why might wind vectors have opposite effect on larval distributions between both lakes? One obvious difference between the two systems studied was the size of captured larvae with the Lake Geneva European whitefish considerably larger compared to Lake Superior cisco larvae (67.8% larger by weight (mean weights = 0.13 g vs. 0.022 g) and 23.7% larger by length (mean lengths = 12.9 mm vs. 15.9 mm). Also, wind speeds were much higher in Lake Superior with wind speeds measured above 47 km.h⁻¹ preceding 7 out of 8 sampling events and maxing out at 67.km hr⁻¹, whereas only one measurement of wind speed was above 20.km hr⁻¹ (20.9 km hr⁻¹) at Lake Geneva. Based on a model predicting larval cisco length by age developed for the south shore of Lake Superior (Oyadomari and Auer, 2008), the larvae we captured inside Grand Portage Bay on 6/14/2018 (13.2 mm, Table 3) were at least 40 days old, coinciding with a predicted hatch date of early May 2018. Larval cisco in Lake Superior can be advected great distances in Lake Superior, with evidence of similar sized larval cisco as the ones captured during this study being advected from spawning grounds in western Lake Superior to the Keweenaw Peninsula (10s of kms away) via currents (Oyadomari and Auer, 2008). Therefore, the southwest winds could have advected larvae into the mouth of Grand Portage Bay prior to the 6/14/2018 sampling event where larval densities inside the bay up until that point had been extremely low. The exact locations of cisco spawning at Grand Portage are unknown (Goodyear et al., 1982) so we cannot say with certainty that cisco spawn and their embryos hatch inside Grand Portage Bay. We can say this seems unlikely given the size of larvae we captured on 6/14/2018 and that few larvae had been captured inside the bay prior to this date.

An alternative explanation about the importance of wind could be that wind impacted our observations of larval densities by decreasing the catch efficiency of our gear. The small larvae might not be able to control their vertical position during high wind speed events in directions that might cause strong wind forcing induced water rotation. For wind events that occurred along the fetch of either lakes, it is possible that larval fish were more vertically dispersed and therefore might not be as present in the surface layer where we were sampling. Or, if the larval fish can hold their position but their prey cannot, they might be following prey to greater depths. However, associations between the vertical distributions of larval coregonids and zooplankton have not been found (Eckmann, 1989; Ventling-Schwank and Meng, 1995; Ylönen et al., 2005). Wind speeds preceding sampling events in Lake Geneva were very low (mean = 2.04 km.hr⁻¹) and were much higher in Lake Superior (mean = 9.17 km.hr⁻¹). Sea conditions on Lake Geneva were calm (<1 m) during our sampling events due to the low wind speeds and were no greater than 0.6 m while sampling on Lake Superior. Larval catches on Lake Geneva were highest on windiest days even when wind direction was along the fetch. Therefore, we feel as if conditions were present where larval fish would have been in the upper surface layer during our sampling. In both the Great Lakes and in European lakes, coregonine larvae are aggregated near the surface during daytime in littoral and pelagic areas (Faber, 1970; Reckahn, 1970; Hoagman, 1973; Viljanen et al., 1995; Selgeby et al., 1978; Hatch and Underhill 1988; Eckmann, 1989; Oyadomari and Auer 2004). Greater than 70% of vendace larvae in Finnish lakes aggregate near the water surface (10–15 cm depth) and 90–96% in the 30 cm layer after hatching (Ylönen et al., 2005; Urpanen et al., 2009). This supports the use of horizontal nets used in our study to investigate larval density (Myers et al., 2008; Karjalainen et al., 2019).

Spring wind patterns might have repercussions on larval cisco survival based on evidence from other fish populations. In Atlantic cod (Gadus morhua), models indicate a strong correlation between wind and age-1 recruitment success with onshore winds having positive effect by keeping eggs and larval fish in nursery areas and offshore winds having a negative effect by advecting eggs and larval fish to less suitable offshore habitat (Churchill et al., 2011). Recruitment success of vendace is reduced during years of heavy wind forcing a month after larvae hatched (Marjomäki, 2004; Marjomäki et al., 2014). Our results suggest that larger coregonine larvae may be more able to maintain position during high wind events and that this generally supports the hypothesized mechanism that Stott et al. (2022) evoked to explain absence and presence of stock structure in Great Lakes cisco populations.

Larval bioenergeticlogs work at three Lake Superior sites (Myers et al., 2014) predicted 10-mm cisco larvae could not grow because of low densities of EZ prey. The coarse level of zooplankton collection (3 stations per site sampled every two weeks) could have missed higher densities of zooplankton as their distributions are patchy. This conclusion of the Myers et al. (2014) study drove us to increase zooplankton sampling at our Lake Superior site to include 10 sites and to increase the frequency of sampling to weekly. Given Lake Superior is ultra-oligotrophic we expected selective pressure to promote maintenance of larval position in food-rich areas, but this was not strongly supported by our data. Densities of EZ at our Lake Superior site were quite low (averages from 1.0 × 10⁻² to 1.7 × 10⁻²ind.L ⁻¹; Tab. 4); well below levels that others have reported finding population-level declines in recruitment. For example, recruitment of Alpine whitefish declined as EZ densities dropped below 10 ind.L⁻¹ in Lake Lucerne, Switzerland (Rellstab et al., 2004). Moreover, recruitment of lake whitefish in the Bay of Quinte, Lake Ontario, declined to very low levels when EZ densities dropped to around 0.4 ind.L⁻¹ (Hoyle et al., 2011). Recent dramatic changes in the stock-recruitment relationship of cisco in western Lake Superior have been hypothesized to result from declining phosphorus levels during the latter half of the 20st century, possibly leading to lower larval survival (Rook et al., 2021). Our estimates of EZ at our Lake Superior site provides further evidence that larval cisco face a daunting challenge of finding food in this ultra-oligotrophic system.

One weakness of our design is we collected larvae in horizontal surface tows, although EZ was measured in standard 25-m vertical tows. Moreover, we did not collect enough paired standard vs. shallow zooplankton tows in Lake Geneva to compare the EZ densities between these habitats. The Lake Superior results did not differ between shallow (5 m to surface) and deeper standard horizontal tows, indicating our measures of EZ in this lake reflected densities available to larvae in surface waters. Depth was a significant driver of larval coregonine catches with higher catches occurring in shallower water however, the shallowest sites were not included until halfway through the study in both lakes. Future efforts to investigate larval densities in these two areas should incorporate more shallow sites, including efforts directly along the shoreline. Additionally, the sampling design for this study was based on Myers et al. (2008) which outlined the best practices for larval sampling in Lake Superior at the time. Due to the hope to compare the two systems, the sampling procedures were mirrored between the two lakes. This, however, meant that sampling protocols used to study larval distributions in European Lakes (Karjalainen et al., 1998; Karjalainen et al., 2002; Viljanen et al., 2002; Urpanen et al., 2009; Leonardsson et al., 2016), were not used which might have led to too broad of a spatial scale of sampling. Sample precision can be increased by taking many small samples or by sampling large volumes of water (Cyr et al., 1992). We sampled large volumes of water per transect (Lake Geneva mean transect volume = 1,354 m⁻³, Lake Superior mean transect volume =160 m⁻³) but did not do replicates of the transects on the same day. Doing a larger number of smaller replicate samples might have increased or improved our ability to understand the fine scale distributions of the larval fish. Follow up research could benefit from more intensive finer scale sampling, at a lake wide scale if feasible, with the possible inclusion of 3D hydrodynamic models similar to Karjalainen et al. (2019). Since this was the first study of larval coregonines in Lake Geneva, we do not have a clear idea of the representativeness of the spawning and nursery habitat at a lake wide scale. For Lake Superior where more research has been done, we are still learning about the spawning behavior of cisco (Shrovnal et al., in review) and we know little about nursery habitat for cisco especially given that our observations do not align with other studies that saw more larval coregonines in greater depths (Oyadomari and Auer, 2008; Myers et al., 2009; 2014; Lucke et al., 2020). Therefore, more research at larger scales would be greatly beneficial to help us understand variation in spawning and nursery habitat across both systems.

We also note that our exploration of factors shaping larval distributions was far from exhaustive. For example, we did not examine how vulnerability to predators across habitats might influence larval distributions (sensu Myers et al., 2014). Although difficult to study, cannibalism in European alpine lakes is suspected to be a major source of mortality of larval coregonines (Eckmann, 2013) as is rainbow smelt (Osmerus mordax) predation on cisco larvae in Lake Superior in some locales (Myers et al., 2009). Coregonine larvae occupation of surface water in the littoral zone may have evolved because their predators (small fishes) are themselves reclusive to avoid being prey. Such a pattern would be consistent with the behavioral trophic cascade hypothesis where the behavior of carnivores can impact vertical distributions of their prey (Romare and Hansson, 2003; Bollens et al., 2011).

Our study of two disparate systems revealed that bathymetric depth is a strong predictor of larval coregonine distributions, suggesting shallow depths were particularly important to the larval life stage in both systems, albeit with the caveat that our study results differ from prior Lake Superior studies showing cisco larvae can occupy great depths. The north-south wind vector and Julian Day were related to larval European whitefish distributions in Lake Geneva and the east-west wind vector was identified as important at our Grand Portage site. The EZ densities we measured in both systems were below levels that others have reported can lead to low larval survival and recruitment failure in other systems. The contrast in the sizes of larvae captured across our study sites, coupled with the relative importance of wind across these two systems, leads us to conclude that larger larvae may be better able to maintain position at preferred locations in the presence of high wind.

Supplementary material

Supplemental Figure 1. Wind rose plots for Lake Geneva (A) and Lake Superior (B) showing wind speed (km.h-1) and direction during the two days prior and first 8 hours of the day of sampling. Mean values of v (north-south wind vector) and u (east-west wind vector) of each sampling event are shown in white box in upper left of each plot. Wind speed measurements of 0 km.hr -1 (49.1% and 0.2% of hourly measurements from Lakes Geneva and Superior, respectively) were omitted. Wind direction and speed data for Lake Geneva was obtained from Sciez weather station (46.331o, 6.3798o), located 10 km southwest of our study area, maintained and operated by © Météo Sciez (http://meteo-sciez.fr). Wind direction and speed data for Lake Superior was obtained from the Rock of Ages (ROAM4) Coastal-Marine Automated Network (C-MAN) station (47.867o, 89.313o), owned and maintained by the National Data Buoy Center of the National Oceanic and Atmospheric Administration (NOAA).

Supplemental Figure 2. Empirical cumulative distribution against fitted distribution functions (Poisson and negative binomial) along with the histograms of larval counts against fitted density functions for Lake Geneva and Lake Superior. These plots were produced using cdfcomp and denscomp from the package ‘fitdistrplus’ and were used to for the most appropriate distribution family to use (Delignette-Muller and Dutang, 2015).

Acknowledgements

The Lake Geneva research was funded with a Visiting Professor Grant awarded by the Université Savoie Mont Blanc, which is one of the supporters of the UMR CARRTEL (INRAE − USMB), whearas the Lake Superior research was made possible by a grant received from the Great Lakes Tribal Initiative − Fishery Management Annex, administered by the U.S. Environmental Protection Agency. This manuscript benefitted by reviews provided by Dr. David Bunnell, Taylor Brown, two anonymous peer reviewers, and assistance from the Guest Editor for the International Journal of Limnology, Dr. Juha Karjalainen. We thank Roger Deschampe Jr., Jim Dahl, Yvette Ibrahim, and Lori Evrard for their help on Lake Superior. Special thanks to INRAE staff including Jean Christophe Hustache and Leslie Laine. This work had support by AnaEE-France and Observatory of Lakes (OLA) boat and technical facilities. Students helped with field and lab processing, often with short notice, including Marine Lemaire, Tristan Flaven, Guillaume Lefebvre, Lisandrina Mari, and Francois Keck. Zooplankton samples were processed by University of Vermont students including Anya Steinhart, Victoria Giacchino, Posy Labombard, and Ben Block. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Appendix A: Development of linear models to predict standard length of zooplankton from total lengths

During quality assurance procedures, we identified incorrect length measurements by one sample processor. The processor had measured lengths to the end of the caudal setae (copepods) and tail spines (cladocerans) for a subset of samples from both lakes. To correct this error, we developed a correction factor to convert total lengths (to end of setae or tail spine) to standard lengths (to end of caudal rami or base of the tail spine). We measured both total and standard lengths from 60–70 individuals from Lake Geneva and 30-40 individuals from Lake Superior, of common zooplankton taxa, from a subset of samples that spanned the entire duration of collections from each lake. We then developed regression models to predict standard lengths from total length measurements and applied these models to the measured (incorrect) lengths of the affected samples to estimate standard lengths. Each model had R² >0.77 (Tab. A1).

Appendix B: Steps taken to define edible zooplankton

B.1 Lake Geneva

Definition of edible zooplankton (EZ) of European whitefish in Lake Geneva was based on the LDS of larvae we caught, understanding of their diets in other alpine lakes in France (Anneville et al., 2007; Anneville et al., 2011) and a reanalysis of zooplankton measured in the stomachs of European whitefish from Lake Annecy, France, from late-February to early-May 2004–2006 (Anneville et al., 2007). The bulk of Lake Geneva larvae captured during our 2016 campaign had LDS scores ≤4 (80%) with 94% ≤19 mm total length. Larval European whitefish exhibit ontogenetic diet shifts (Anneville et al., 2007), and reanalysis of their data showed LDS 1–2 consumed mainly cyclopoid copepod copepodites and adults <0.7 mm with some use of nauplii. The LDS 3–4 also consumed cyclopoid copepods and adults < 0.7 mm, and ingested cladocerans up to 0.6 mm in length. Thus, EZ of Lake Geneva larvae were defined as cyclopoid copepods ≤0.7 mm, cladocerans ≤0.6 mm, and nauplii. Calanoid copepods, regardless of size, were excluded (see Tab. B1 for comparison between EZ vs all zooplankton length ranges and density estimates).

B.1 Lake Superior

We determined the zooplankton length cut-off for edible prey on Lake Superior larvae based on our genetic assignments (Table 4), our length measurements, results of a recent Lake Superior larval diet study (Lucke et al., 2020) and estimates of gape width of coregonine larvae of varying lengths (Davis and Todd, 1998). Gape width of larval fish effects the upper limits of consumed prey sizes, although most prey consumed are generally ~0.2 mm smaller than the maximum across the range of larval lengths (see Fig. 4 in Schael et al., 1991).

Larval coregonines collected between 5/14/2018 and 6/13/2018 in the Apostle Islands of Lake Superior's south shore consumed invertebrate eggs, copepod nauplii, and calanoid copepod copepodites and adults, but not cyclopoid copepods and cladocerans (Lucke et al., 2020). Lachance et al. (2021) later assigned larvae from the Lucke et al. (2020) study to species and showed most collected prior to 6/11/2018 were cisco (see their Fig. 3). Because the larvae we caught at our Lake Superior site were also primarily cisco (97%) and 88% of all larvae we captured were ≤14 mm, we opted to limit EZ to individuals <0.75 mm. This was based on the work of Davis and Todd (1998) showing 14-mm cisco larvae have a mean gape of 0.94 mm (see their Tab. 2). The 0.75 mm maximum prey size is a recognition that the bulk of zooplankton consumed by larvae in the study of Schael et al. (2011) were < ~0.2 mm smaller than the maximum based on gape (see their Fig. 3). To recap, EZ for Lake Superior larvae were defined as copepod nauplii, calanoid copepodites and adults <0.75 mm; cladocerans and cyclopoid copepods, regardless of size, were excluded. (see Tab. B1 for comparison between EZ vs all zooplankton length ranges and density estimates).

References

All Tables

All Figures

Fig. 1 Locations of zooplankton (black circles) and larval tow stations (black lines) sampled during our Lakes Geneva (2016) and Superior (2018) campaigns. The Lake Geneva temperature logger at the laboratory port break wall is labeled INRAE (the National Research Institute for Agriculture, Food and the Environment laboratory), and the surface mini-logger station on Lake Superior is labeled Mini-logger.

In the text

Fig. 2 Surface temperatures (Celsius, °C) versus calendar date (MM/DD) at Lake Geneva sites sampled during 2016 and Lake Superior sites during 2018. The Lake Geneva temperatures are the averages from 0–1 m depths from water column profiles and the thermistor positioned at 1m depth at the National Research Institute for Agriculture, Food and the Environment laboratory (INRAE) laboratory break wall (Fig. 1). The Lake Superior temperatures are the average of 24-hourly recordings measured by a data logger deployed at the surface at a single site (Mini-logger; Fig. 1).

In the text

Fig. 3 Larval coregonid natural log-transformed densities (ind.103m−3) and total lengths (mm) for each larval fish by sampling date in Lake Geneva and Lake Superior. Shallow site (<10 m bathymetric depths) larval densities (Lake Geneva site 0 and Lake Superior sites B1, B2 and B3, Fig. 1) and lengths are shown with open diamonds. Before transforming the data, half of the smallest non-zero value were added for each lake (0.3 for Lake Geneva and 1.6 for Lake Superior) due to the presence of zeros (Leith et al., 2010) which places zero values at −1.20 for Lake Geneva and 0.47 for Lake Superior after transformation. The black line within the box marks the median, the boundary of the box closest to zero indicates the 25th percentile, and the boundary of the box farthest from zero indicates the 75th percentile. Whiskers represent 1.5 * IQR from the bottom and top of the boxes, where IQR is the inter-quartile range or distance between the first and third quartiles. Data beyond the whiskers are outliers.

In the text

Fig. 4 Edible zooplankton densities (ind .L−1) for larvae based on 64-μm mesh tows at the Lake Geneva stations as a function of bathymetric depth range and sampling date. Stations were sampled between 4/4/2016 and 4/21/2016. Stations were binned as follows: <10 m bathymetric depths = site 0; 20–55 m bathymetric depths = sites 1–4; 60–100 m bathymetric depths = sites 5–8 (see Fig. 1). Stations <10 m were added on 4/14/2016. The 4/4/16 samples were collected using a closeable net instead of a double Bongo net, so density estimates were adjusted based on a conversion factor we developed (see Results section). Edible zooplankton on Lake Geneva was defined as cyclopoid copepods ≤0.7 mm, cladocerans ≤0.6 mm, and nauplii (see Appendix B for details).

In the text

Fig. 5 Edible zooplankton densities (ind L−1) for larvae based on 64-μm mesh tows at the Lake Superior stations by bathymetric depth range and sampling date. Stations were sampled between 5/3/2018 and 6/20/2018. Stations were binned as follows: < 10 m bathymetric depths = B1; 20-55 m bathymetric depths = sites S1-S2 and M3-M6; 60-100 m bathymetric depths = sites D7-D10 (see Fig. 1). The < 10 m zooplankton station was added on 5/31/2018. Edible zooplankton on Lake Superior was defined as copepod nauplii, calanoid copepodites and adults <0.75 mm; cladocerans and cyclopoid copepods, regardless of size, were excluded (see Appendix B for details).

In the text

Fig. 6 Generalized additive model predicting surface larval counts in Lake Geneva (left panel) and larval counts with volume as an offset in Lake Superior (right panel) as a function of explanatory variables. Circles are larval density observations and shaded areas are 95% confidence intervals. Jittering was applied when plotting due to overlapping data points. A break in the y-axis was applied for Lake Superior due to a high larval density observation and horizontal lines are provided to accentuate the break.

In the text

	Fig. 7 Observations for predictor variables not included in the top models for each lake vs. larval density. Jittering was applied when plotting due to overlapping data points.
In the text