Reconstructing half a century of coregonine recruitment reveals species-specific dynamics and synchrony across the Laurentian Great Lakes

Abstract

Understanding how multiple species and populations vary in their recruitment dynamics can elucidate the processes driving recruitment across space and time. Lake Whitefish (Coregonus clupeaformis) and Cisco (C. artedi) are socioecologically important fishes across their range; however, many Laurentian Great Lakes populations have experienced declining, poor, or sporadic recruitment in recent decades. We integrated catch and age data from 38 long-term surveys across each of the Great Lakes and Lake Simcoe, resulting in a combined time series spanning 1960–2019. We estimated Lake Whitefish and Cisco year-class strength (YCS) in each lake using longitudinal mixed-effects regressions of relative cohort abundance. We subsequently quantified interspecific, spatial, and temporal synchrony in YCS using correlation and dynamic factor analyses. Lake Whitefish YCS was positively spatially synchronous on average, and YCS in all six lakes was elevated during the 1980s–1990s. In contrast, Cisco YCS was sporadic, not spatially synchronous, and highly variable around long-term, lake-specific means. YCS was not synchronous between species in any lake. Collectively, our analyses demonstrate that these species exhibit differential recruitment dynamics that may be regulated by species-specific factors. Results from this study can be leveraged in future research on the causes and consequences of cross-species, cross-basin recruitment variability.

year-class strength, time series analysis, Coregonus, Lake Whitefish, Cisco, Simcoe

Introduction

Recruitment is a fundamental biological process governing the status and trajectory of populations. Fish recruitment—defined as the number of individuals from a given cohort that survive to a focal stage (e.g. sexual maturity, reaching a size vulnerable to fishing gear; Krabbenhoft et al. 2023)—is often thought to be regulated by processes acting on early life stages (Houde 2008). Understanding recruitment has proven to be one of the most enduring challenges in fisheries science because recruitment is inherently variable and can be regulated by many interacting factors, the relative importance of which can vary spatially and temporally (Ricker 1954, Myers 1998, Subbey et al. 2014, Szuwalski et al. 2015, Munch et al. 2018). Recruitment drivers for a given species may differ among populations due to variation in local ecosystem context and spawning stock sizes (Brosset et al. 2019, Krabbenhoft et al. 2023) and may be nonstationary within populations due to environmental change and shifting demographics (Walters 1987, Bunnell et al. 2006, Vert-pre et al. 2013, Feiner et al. 2015). Because of this complexity, many predictive relationships of recruitment based on single populations perform poorly when applied outside of the original sample frame (Myers 1998).

Understanding the extent to which recruitment dynamics (e.g. frequency of strong year-classes) vary across populations spanning biophysical gradients can enhance our ability to predict future dynamics and responses to ecosystem change. For example, determining the degree of spatial synchrony can inform how populations will respond to continued climatic change. Climatic processes operating across large spatial scales can spatially synchronize recruitment across discrete populations (Moran 1953, Koenig 2002), whereas populations may be asynchronous if local processes are more influential (Rogers and Schindler 2008). Further, recruitment periodicity can be regulated by intrinsic and extrinsic factors, such as life history traits (Winemiller 2005), spawning stock biomass (Ricker 1954, Beverton and Holt 1957), environmental variability (Vert-pre et al. 2013), and fishing pressure (Anderson et al. 2008). Investigating recruitment patterns (e.g. consistent interannual recruitment, strong periodicity, decadal shifts) can help set expectations for stock productivity under differing contexts and is particularly pertinent for predicting future returns on fished populations (Kuparinen et al. 2014) or extirpation risk for declining populations (Durham and Wilde 2008). Management strategies are most likely to be effective when tailored to expectations for a given species, but even species-specific expectations may need to change through time due to nonstationarity (Feiner et al. 2015, Marcek et al. 2021). Improving understanding of the causes and consequences of recruitment variability is critical for fisheries managers and stewards to develop effective strategies for sustaining fisheries in an era of accelerating ecosystem change (DeVanna Fussell et al. 2016).

The Laurentian Great Lakes of North America (hereafter Great Lakes) are an ideal setting for investigating cross-system recruitment dynamics. The immense spatial scale of the Great Lakes results in physical regimes analogous to those of inland seas (Ludsin et al. 2014, Sterner et al. 2017). Consequently, the recruitment dynamics of many Great Lakes fishes operate more similarly to that of fishes in marine systems (Janssen et al. 2014, Ludsin et al. 2014, Pritt et al. 2014), where broad-scale climatic and physical processes regulating early life stage survival act to drive recruitment and synchronize dynamics across populations (Myers et al. 1997), than recruitment dynamics of fishes in smaller freshwater lakes, where local factors that influence juvenile survival can be more important (Houde 1994). For example, Great Lakes studies have documented recruitment synchrony both among populations of the same species with limited dispersal (Bunnell et al. 2010, Honsey et al. 2016, Warren et al. 2024) and among species with diverse life histories (Bunnell et al. 2017), highlighting the importance of broad-scale processes for regulating recruitment. Further, each of the five Great Lakes—Superior, Michigan, Huron, Erie, and Ontario—is uniquely situated across gradients of physical regimes, biological community structures, and anthropogenic influences (Ives et al. 2019). The lakes share many fundamental characteristics related to climate, geology, and biogeography, but biophysical conditions can substantially differ across lakes, including latitudinal variation in weather, primary productivity, habitat degradation, the prevalence of non-native species and their subsequent effects on food-web structure, and fishing intensity (Allan et al. 2013, Dove and Chapra 2015, Ives et al. 2019). Comparing recruitment dynamics across these systems can help to clarify context-dependent population dynamics and responses to ecosystem change.

Coregonines (Salmonidae Coregoninae) are ecologically, economically, and culturally important coldwater fishes distributed across the northern temperate zones of Eurasia and North America (Anneville et al. 2015). Two coregonine species, Lake Whitefish (Atikameg in Anishinaabemowin; Coregonus clupeaformis) and Cisco (Otoonapii; C. artedi), are central to ongoing efforts to sustainably manage, conserve, and restore native fish populations across the Great Lakes region (Bunnell et al. 2023). Lake Whitefish and Cisco contribute to native fish diversity (Koelz 1929, Eshenroder et al. 2016), support food web function (Stockwell et al. 2014, Matthias et al. 2021), are culturally significant (Gobin et al. 2022, Duncan et al. 2023), and are commercially harvested in the Great Lakes (Great Lakes Fishery Commission 2022). Contemporary populations of both species are in varied states of recovery following fishery collapses throughout the late 1800s to mid-1900s and subsequent decades of harvest regulation, habitat restoration, and invasive Sea Lamprey (Petromyzon marinus) control—except in Lake Erie, where Cisco are considered extirpated (Eshenroder et al. 2016). Disentangling the processes driving recruitment is critical for effective management and stewardship of Lake Whitefish and Cisco populations (Zimmerman and Krueger 2009, Council of Lake Committees 2018, Ebener et al. 2021). Notably, it is unclear why recruitment is sporadic or declining for many populations but consistent or improving for others. For example, a few Lake Whitefish populations (e.g. Ransom et al. 2021, Carl et al. 2024) appear to be exceptions to widespread declines in recruitment since the early 2000 s (Ebener et al. 2021), and recruitment trends for some co-occurring Lake Whitefish and Cisco populations have recently diverged (e.g. Claramunt et al. 2019, Brown et al. 2022, Modeling Subcommittee 2022). These species share generally similar early life histories (Brown et al. 2023) that differentiate after the larval stage when Lake Whitefish transition to demersal habitats (Reckahn 1970) and Cisco remain pelagic (George 2019), which could result in species-specific recruitment bottlenecks in post-larval life stages. Increased understanding of the processes regulating the dynamics of these two species could help clarify the causes of declining, sporadic, and divergent recruitment and can inform management (Ebener et al. 2021, Bunnell et al. 2023).

Cross-lake comparisons can elucidate underlying patterns and important differences in recruitment dynamics for Lake Whitefish and Cisco. For one, differences in biophysical conditions among lakes could be driving variable recruitment success across the Great Lakes (Brown et al. 2024). Previous studies suggest that Cisco exhibit spatial synchrony in recruitment within and among lakes (Rook et al. 2012, Myers et al. 2015, Weidel et al. 2021) whereas Lake Whitefish recruitment is asynchronous within lakes (Zischke et al. 2017). However, cross-lake studies of recruitment to date have generally been limited to the Upper Great Lakes (Superior, Michigan, and Huron), which are more similar to each other (Barbiero et al. 2012) than to the Lower Great Lakes (Erie and Ontario), potentially masking important differences in dynamics exhibited across a wider range of biophysical gradients. Previous studies of Lake Whitefish and Cisco synchrony have also applied different methods at varying spatial scales, complicating direct inferences among studies. More broadly, comparisons between the Great Lakes and smaller inland lakes could help resolve the scale at which factors that regulate recruitment operate (e.g. local, lake-wide, regional). Recruitment dynamics can also differ between lakes based on the scale of the lakes themselves (e.g. Krabbenhoft et al. 2023), as local factors are likely more important for small lakes than for large lakes (Rudstam 1984, Myers et al. 1997, Janssen et al. 2014). Lake Simcoe (Ontario, Canada) is geographically proximate to the Great Lakes and supports fisheries for Lake Whitefish and Cisco (MacCrimmon and Skobe 1970, Dunlop et al. 2019). The Lake Simcoe ecosystem has undergone similar environmental changes to the Great Lakes (e.g. Goto et al. 2020) but is more productive (North et al. 2013) and relatively small in surface area (744 km²) compared to the Great Lakes (smallest surface area = 19 000 km² for Lake Ontario). For these reasons, Lake Simcoe provides an opportunity to further investigate how recruitment dynamics converge or diverge across lakes that vary in scale and biophysical conditions.

The overarching goal of this study was to characterize and compare Lake Whitefish and Cisco recruitment dynamics among six lakes: Superior, Michigan, Huron, Erie, Ontario, and Simcoe. To accomplish this, we generated standardized year-class strength (YCS) indices for each species and lake by applying longitudinal mixed-effects regressions of relative cohort abundance data from 38 long-term surveys with catch and age data. As in He et al. (2023), we define YCS as the relative abundance of a fish year-class as derived from repeated measures across multiple ages and sampling years; importantly, this definition differentiates YCS from recruitment indices based on a given age class as observed in a single year. We also evaluated the utility of two related metrics of relative abundance in reconstructing YCS: catch-at-age and age-structured catch-per-unit-effort (CPUE). We subsequently quantified synchrony in YCS between species, among lakes, and through time using correlation and dynamic factor analyses. We hypothesized that: (1) Lake Whitefish and Cisco recruitment is primarily regulated by processes acting on early life stages and, given similarities in their early life histories, the two species share important drivers of recruitment that act to synchronize interspecific YCS; (2) climatic processes are important drivers of early life stage survival and act to synchronize YCS across the Great Lakes region, but differences in biophysical conditions among lakes determine the degree of cross-lake synchrony; and (3) concurrent ecosystem changes across the Great Lakes have induced common temporal trends in YCS among lakes for each species. Results from this study advance our understanding of the recruitment dynamics of Lake Whitefish and Cisco over half a century, inform ongoing management and stewardship frameworks, and can be leveraged in future research on the causes and consequences of cross-species, cross-basin recruitment variability.

Materials and Methods

Studies of recruitment across populations and systems must be robust to differences in sampling design among surveys and through time. Because the Great Lakes are large systems with many management agencies and jurisdictions, there is a wealth of information on Lake Whitefish and Cisco catches, but surveys can vary considerably in design. Fishery-dependent surveys (FDS) tend to be prioritized for fisheries monitoring and observe large numbers of fish. However, FDS often employ nonrandom sampling designs and are highly size-selective; further, harvest rates do not necessarily reflect overall trends in populations (Hilborn and Walters 1992). In contrast, fishery-independent surveys (FIS) typically employ standardized sampling designs that seek to provide catch rates and biological information in proportion to the overall population. However, many FIS often span relatively short time periods and can have low sample sizes in some years, especially when population abundance is low. Few FIS in the Great Lakes include coregonines as target species and therefore may not prioritize full biological workup of captured Lake Whitefish and Cisco. Ultimately, methods that integrate multiple FIS and FDS with varying sampling designs within a single modeling framework can provide robust estimates of recruitment (Pennino et al. 2016, He et al. 2023).

Data compilation and harmonization

We compiled long-term catch and age data describing Lake Whitefish and Cisco populations within each of the Great Lakes and Lake Simcoe (Appendix Table A1). These data originated from 24 FIS and 14 FDS conducted by state, provincial, Tribal, and United States federal fisheries agencies across the region (Supplementary Material S1). Our goal was to estimate YCS for individual cohorts as far into the past as possible; therefore, we selected surveys that targeted and/or frequently encountered adult or juvenile Lake Whitefish or Cisco, had robust age data spanning long time periods, and had appropriate sampling designs for our inferential and analytical goals (e.g. representative fishing gears). Prospective surveys were identified through the Coregonine Restoration Framework (Bunnell et al. 2023) and by leveraging the authors’ knowledge of available data. For each survey, we compiled data on sampling effort, total catch, and biological data for individual Lake Whitefish and Cisco. We also compiled fishery-specific harvest statistics associated with each FDS. Though available for some management units (e.g. Modeling Subcommittee 2022), we chose not to use estimated numbers of recruits derived from stock assessment models as a data source to avoid structural assumptions (e.g. stock-recruitment relationships; Ricker 1954) that may or may not be accurate (Brooks and Deroba 2015).

We accounted for differing size selectivity and calculated standardized metrics of sampling effort among sampling gears. Because age and size are often correlated in fishes, the choice of sampling gear can influence the observed age composition (Guy and Brown 2007). Single surveys commonly encompassed multiple sampling gears (e.g. commercial monitoring from both gill nets and trap nets), so we treated each sampling gear separately for each survey to account for differing size selectivity. Surveys included in our analysis were typically conducted using gill nets and/or impoundment nets (i.e. trap nets, pound nets, and hoop nets) with a few exceptions: two recreational creel surveys and one fishery-independent electrofishing survey. We further categorized gill net data by mesh size to account for differing size selectivity: “large” mesh (≥ 4 in or 101 mm stretched mesh), “small” mesh (< 4 in), or “graded” mesh (i.e. multi-panel meshes designed to select for a wide range of sizes). We standardized sampling effort for FIS as one km of gill net or one impoundment net lift. We did not use available metrics of sampling effort from FDS to avoid assumptions of using fishing effort-derived catch-per-unit-effort (CPUE), as they can be unreliable metrics of abundance over time (Harley et al. 2001). Similarly, we did not have reliable metrics of effort from the fishery-independent electrofishing survey and thus only used catch information. Our standardized metrics of effort did not account for sampling duration, as we expected soak time to have a non-linear and density-dependent effect on catch rates (e.g. Hansen et al. 1998), and duration was not commonly recorded among surveys. Lastly, gill nets have shifted from multifilament to monofilament material through time (e.g. Rook et al. 2022b), and catchability can differ between the two net materials, particularly with increasing water clarity (Smith et al. 2022). While we did not explicitly account for this change, our model structure was designed to capture changes in catchability through time (details below).

We accounted for potential measurement error in fish age data, which were interpreted from otoliths, scales, or fin rays by the original data collectors. Otoliths and fin rays have recently become common for interpreting ages because scales underestimate ages of Lake Whitefish (Muir et al. 2008a, 2008b) and Cisco (Rudstam 1984, Yule et al. 2008), with measurement error increasing with fish age, though each choice of age structure and preparation technique has inherent measurement error (McKeefry et al. 2023). However, scales continue to be used today, especially for smaller-bodied individuals, and historical age information is almost universally based on scales. We accounted for potential age underestimation for older fish by limiting the maximum ages derived from scales included in our analyses. Lake Whitefish scale ages begin significantly deviating from otolith ages between ages 5–8 years, depending on the population (Muir et al. 2008a, 2008b, Herbst and Marsden 2011). Scale ages from Lake Superior Cisco begin deviating from otolith ages starting at age 5 (Yule et al. 2008), but information from other populations is limited. However, there is no single threshold that will eliminate aging bias because the discrepancy among hard structures is cumulative with fish age. Many surveys begin observing fish at older ages (e.g. age 6) due to gears that select for larger-bodied fish; therefore, the choice of maximum age to include for scale-aged individuals can lead to large losses of available historical and contemporary data. We chose age 7 as a maximum age for scale-aged individuals for both species to balance the trade-off between confidence in aging accuracy and data availability. Importantly, this choice had minimal impacts on our results (Supplementary Material S2, Supplementary Figs. S1, S2). When age structures were not recorded, we conservatively assumed scales were used.

For each survey, we summarized species-specific catch (expressed in numbers of fish) and age data to estimate annual catch-at-age and age-structured CPUE. We estimated both catch-at-age and CPUE for FIS with standardized effort information, whereas we only estimated annual catch-at-age for surveys without reliable effort information (i.e. all FDS, electrofishing survey). We estimated annual catch-at-age for each survey, including observed but unaged fish, using observed age compositions and extrapolating to total catch each year. We also calculated annual CPUE by dividing catch-at-age by the total standardized sampling effort each year. For FDS, we used harvest statistics to calculate catch-at-age because FDS biological monitoring often uses defined subsampling quotas (e.g. up to 50 fish aged each year, regardless of harvest). We converted fishery-specific harvest reported in total biomass to estimated number of individuals using the mean weight of sampled fish. We converted dressed weights to round weights when necessary using standard conversion factors for each species (Great Lakes Fishery Commission 2022). When possible, we estimated the number of individuals harvested based on weight data specific to each survey, year, and sampling gear to maximize the accuracy of our estimates, but biological data were not available for all combinations. When weight data were unavailable, we borrowed weight estimates from (a) each year for a given FDS across all fishery-dependent sampling gears (19% of sampling years across all surveys and sampling gears), (b) all FDS biological data within each lake for a given year (4%), or (c) all FDS biological data for a given year (<1%), in order of preference. As with FIS, we estimated the total annual catch of each age class using known age composition and extrapolating to the total harvest each year. We used observed catch numbers directly when total harvest data were unavailable (7%, not including one creel survey for which we exclusively used observed catch due to a lack of reliable harvest data).

YCS estimation

We estimated YCS using longitudinal linear mixed models of catch-at-age and age-structured CPUE. This approach uses the observed catch of multiple age classes through time to estimate the relative abundance of each cohort (He et al. 2023). We chose this approach because it provides several advantages for quantifying trends in recruitment across lakes and species. Specifically, this method incorporates repeated observations of multiple age classes through time to estimate the YCS of each cohort and should therefore be more robust to sampling bias than single-age recruitment indices. The method is also flexible and allows for multiple sources of information to be combined into a single model, including both FIS and FDS without available indices of effort. This method can also be used with only the subset of age classes that are reliably caught and accurately aged, helping to mitigate issues associated with old, rare individuals and aging bias incurred by some approaches (e.g. Honsey et al. 2016, Zischke et al. 2017). Lastly, this method avoids assumptions inherent to stock-recruitment relationships. We used all available catch and age data in each lake to estimate lake-wide YCS for each species. This lake-wide approach allowed us to comprehensively reconstruct Lake Whitefish and Cisco YCS over multiple decades in each lake.

Our YCS model structure was adapted from those developed by Honsey et al. (2020) and He et al. (2023) to recognize the complex reality that different surveys may observe the relative abundance of each year-class at different age ranges and sampling years. To estimate YCS, we applied the following regression model separately to our two metrics of relative abundance, catch-at-age and age-structured CPUE:

$$\begin{eqnarray} \ln \left( {\it{Inde}{{x}_{i,j,k,m}}} \right) = {{{\rm{\beta }}}_i} + {{a}_{j,m}} + {{c}_{k,m}} + {{{\rm{\varepsilon }}}_{i,j,k,m}}, \end{eqnarray}$$

(1)

where Index represents either catch-at-age or age-structured CPUE of year-class i at age j from sampling year k and data source m; |${{{\rm{\beta }}}_i}$| is a fixed-effect parameter for year-class (i.e. YCS); |${{a}_{j,m}}$| is a normally distributed random intercept among each age j within each data source m; |${{c}_{k,m}}$| is a normally distributed random intercept among each sampling year k within each data source m; and |${{{\rm{\varepsilon }}}_{i,j,k,m}}$| is the residual error following a normal distribution. Each data source is a single combination of survey and sampling gear. YCS is estimated as the coefficient on each factor level (i.e. cohort) of |${{{\rm{\beta }}}_i}$|⁠. The |${{a}_{j,m}}$| random effect accounts for differing size-selectivity across ages, surveys, and sampling gears. The random effect |${{c}_{k,m}}$| accounts for changes in catchability, sampling design, and sampling effort through time, thereby allowing catch to be modeled directly without specifying sampling effort (He et al. 2023). While this approach has separately been applied to independent datasets of catch-at-age (He et al. 2023) and CPUE (Honsey et al. 2020), the modeling framework has yet to be implemented in parallel. Here, we do so to evaluate the utility of the two metrics of relative abundance in estimating YCS. We z-score normalized (i.e. to |$N( {0,1} )$|⁠) our input response data (i.e. |$\ln ( {\it{Inde}{{x}_{i,j,k,m}}} )$|⁠) within each survey prior to modeling. We did this to equally weight the variability in catch across surveys, as annual catch can vary in magnitude based on the nature of each survey (e.g. targeted commercial harvest monitoring programs versus fish community surveys). Without this standardization, YCS estimates for each lake and species would not be directly comparable because the estimated coefficients are dependent on the range of observed data, and differing abundances affect the number of recruits from each cohort. All statistical analyses were conducted in R version 4.3 (R Core Team 2023). We fit separate models for each combination of species and lake using the lme4 package (Bates et al. 2015). Model fits were validated through visual inspection of the residuals on both fixed and random components of all models using the performance package (Lüdecke et al. 2021).

Prior to analysis, we refined our dataset to ensure all data informing the models were reliable sources of YCS information through time. Following He et al. (2023), we restricted ages to those that were fully selected to each gear type for each species in each lake based on observed age distributions over time (Appendix Table A2). Inclusion of younger individuals not fully vulnerable to a gear can underestimate YCS, whereas the inclusion of older individuals can introduce confounding effects of post-recruitment mortality. We therefore limited the analysis frame for each cohort to a three-year window, starting with the first age fully selected to each gear. We used a shorter, two-year window instead for the few instances where strong size-selectivity resulted in only two ages being consistently caught (Table A2). We used a global maximum age of 8 to limit the influence of older age classes and excluded cohorts that were observed in only one sampling year. Our models therefore followed the abundance of a given cohort through time within the range of age 3–8, depending on the lake and sampling gear. We evaluated the degree to which our choice of age ranges impacted our results, both in general (Supplementary Figs. S3-S4) and in the context of shifts in Lake Whitefish size-at-age over time (Supplementary Fig. S5; Mohr and Nalepa 2005). We accounted for this by ensuring our age ranges were fully selected to the gears both before and after observed growth declines based on sampled age distributions for each lake. Importantly, our YCS estimates were robust to the choice of age ranges (Supplementary Material S2, Supplementary Figs. S1-S5).

Interspecific, spatial, and temporal synchrony

We used our standardized YCS estimates to investigate multiple aspects of synchrony between species, among lakes, and through time. First, we examined whether YCS was synchronous between Lake Whitefish and Cisco within each lake. Second, we quantified spatial synchrony in YCS across the Great Lakes and Lake Simcoe for each species. Third, we identified common trends in YCS through time among lakes for each species using dynamic factor analysis (DFA). We used the YCS estimates derived from catch-at-age data for all synchrony analyses, as these models included all available data, allowed for the analysis of longer time series, and produced similar YCS estimates to the CPUE model (see section “Results”).

We estimated interspecific synchrony between Lake Whitefish and Cisco using correlation analysis. We calculated Pearson’s product-moment correlation coefficients (r) between Lake Whitefish and Cisco YCS in each lake. Negative r values indicate negative synchrony (i.e. opposing trends), values near 0 indicate asynchrony, and positive r values indicate positive synchrony (i.e. agreeing trends). We evaluated evidence of interspecific synchrony (i.e. whether r significantly differed from 0) in each lake by conducting two-sided, one-sample t-tests at an α = 0.05 significance level.

We also used correlation analysis to estimate spatial (i.e. among-lake) synchrony for each species. We calculated r in YCS among all pairwise combinations of lakes and averaged all pairwise r values for each species to obtain the mean correlation in YCS across lakes (⁠|$\bar{r}$|⁠). We evaluated evidence of synchrony for each species across all lakes (i.e. whether |$\bar{r}$| significantly differed from 0) by conducting a two-sided, one-sample t-test at α = 0.05. We also tested for significant correlations between all pairwise combinations of lakes. Because we conducted 15 and 10 simultaneous comparisons for Lake Whitefish and Cisco, respectively, we corrected for multiple hypothesis testing using the false discovery rate (Benjamini and Yekutieli 2001) and deemed the adjusted P-values significant at α = 0.05.

Lastly, we used DFA to identify common temporal trends among the lake-specific YCS time series. Briefly, DFA is an autoregressive state-space model wherein each time series is treated as a single observation process of one or more underlying latent trends, all of which can be linearly combined to reconstruct the observations (Zuur et al. 2003). The overarching goal of DFA is dimension reduction, where the fewest number of latent trends are identified to minimize information loss. We implemented DFA independently for each species using the MARSS package (Holmes et al. 2012), which models each trend as a random walk. Prior to analysis, we z-score normalized each lake-specific YCS time series and evaluated whether autoregressive model assumptions were met (e.g. by estimating autocorrelation functions). We used model selection to identify the best supported number of latent trends and error structure among the time series. Potential error structures included the same variance and covariance, same variance with no covariance, differing variance with no covariance, and differing variance and covariance (Zuur et al. 2003). Candidate models that failed to converge after 1 000 000 iterations were removed from further analysis. We used Akaike's information criteria corrected for small sample sizes (AICc) to assess the relative support of fitted DFA models (Akaike 1974) and deemed models with ΔAICc ≤ 2 to be most highly supported (Burnham and Anderson 2002). We validated the best supported model for each species by visually inspecting model residuals before extracting the estimated latent trends, factor loadings, and predicted model fits for interpretation.

Results

We estimated YCS for Lake Whitefish and Cisco cohorts spanning 1956–2015 across the Great Lakes and Lake Simcoe (Figs 1 and 2). Our analyses were informed by 38 long-term surveys with a mean (± 1 standard deviation) time series duration of 23.6 ± 12.6 years (Table A1). Lake Whitefish YCS was estimated for lakes Superior (n = 11 surveys), Michigan (n = 10), Huron (n = 8), Erie (n = 2), Ontario (n = 2), and Simcoe (n = 3). Cisco YCS was estimated for lakes Superior (n = 9 surveys), Michigan (n = 3), Huron (n = 9), Ontario (n = 2), and Simcoe (n = 3). We were unable to include Cisco from Lake Erie in our analysis because they are considered extirpated (Eshenroder et al. 2016).

Figure 1.

Year-class strength (YCS) time series for Lake Whitefish in the Laurentian Great Lakes and Lake Simcoe. YCS was estimated independently for each study lake (columns) using two different metrics of relative abundance: catch-at-age and age-structured catch-per-unit-effort (CPUE). Shaded areas depict 95% confidence intervals. Y-axis units are in standard deviations, where 0 (dashed line) denotes mean YCS relative to each respective model.

Figure 2.

Year-class strength (YCS) time series for Cisco in the Laurentian Great Lakes and Lake Simcoe. YCS was estimated independently for each study lake (columns) using two different metrics of relative abundance: catch-at-age and age-structured catch-per-unit-effort (CPUE). Shaded areas depict 95% confidence intervals. Y-axis units are in standard deviations, where 0 (dashed line) denotes mean YCS relative to each respective model. Cisco are considered extirpated from Lake Erie and therefore are not included in this analysis.

We did not identify structural problems in any fitted models, nor did we detect any major problems across global maximum ages (Supplementary Figs. S1–S2), maximum scale ages (Supplementary Figs. S3–S4), or changes in size-at-age through time (Supplementary Fig. S5). Overall, the catch-at-age and age-structured CPUE models produced similar YCS estimates for most of the time series (Figs. 1 and  2). Advantageously, the models based on catch-at-age (FIS and FDS data) were able to estimate YCS for 13 more cohorts than those based on CPUE (FIS only). Most estimates from the two response variables agreed in direction and magnitude, thus providing equivalent interpretations for YCS variability through time. However, some portions of the modeled time series deviated from one another, resulting in conflicting interpretations of YCS (e.g. strong versus weak) between response variables. To investigate the factors contributing to these discrepancies, we used two approaches to evaluate the utility of the two metrics of relative abundance: (1) separately modeling catch-at-age and CPUE using the same FIS data with known indices of effort, and (2) modeling catch-at-age using only FIS or FDS data, respectively. Many of the discrepancies between the two metrics of relative abundance were resolved using the first approach (Supplementary Fig. S6), indicating that the model’s random effect was successfully capturing variation in sampling effort through time. However, when we compared YCS estimates from FIS versus FDS data alone, many of the original discrepancies remained or became magnified (Supplementary Fig. S7). We concluded that the discrepancies between the two response variables were primarily attributable to occasional differences in YCS trends between FIS and FDS data sources and, to a lesser degree, uncaptured variation in sampling effort. Consequently, we proceeded with the more robust catch-at-age model results for all downstream analyses.

Lake Whitefish and Cisco cohorts commonly exhibited differential YCS in each lake and we did not detect significant interspecific synchrony in any of the six lakes (all |$| r |$| ≤ 0.29, P > 0.05; Fig. 3). While estimated YCS was similar for some cohorts of both species (e.g. 2004, 2008, and 2012 were strong year-classes for both species in Lake Simcoe), agreement between species was relatively rare across the study period. Across all lakes, Cisco YCS was highly variable around lake-specific long-term means, and the frequency of strong year-classes was sporadic. In contrast, Lake Whitefish often exhibited periods of relatively consistent YCS within each lake.

Figure 3.

Year-class strength (YCS) time series and interspecific synchrony analyses for Lake Whitefish and Cisco within each study lake (panels). Inset boxes report the estimated correlations (r), degrees of freedom (df), and P-values for associated hypothesis test (see text for details). YCS was estimated using catch-at-age. Shaded areas depict 95% confidence intervals. Y-axis units are in standard deviations, where 0 (dashed line) denotes mean YCS relative to each model. Cisco are considered extirpated from Lake Erie and therefore are not included in this analysis.

On average, Lake Whitefish YCS was positively synchronous among lakes (⁠|$\bar{r}$| = 0.29, P < 0.001), whereas Cisco YCS was asynchronous among lakes (⁠|$\bar{r}$| = 0.02, P = 0.717). All significant pairwise correlations between lakes were positive for Lake Whitefish (Table 1). Estimated spatial synchrony was strongest for Lake Whitefish between lakes Huron and Michigan (r = 0.77, P < 0.001). We also detected positive synchrony in Lake Whitefish YCS between Lake Simcoe and the Great Lakes. In contrast to Lake Whitefish, no pairwise correlations of Cisco YCS were significantly different from 0 (all |$| r |$| ≤ 0.49, P ≥ 0.452), suggesting that Cisco YCS was asynchronous among lakes.

Table 1.

Spatial synchrony analysis using pairwise Pearson correlation coefficients (below diagonal) and corresponding P-values adjusted for multiple comparisons (above diagonal).

	Superior	Michigan	Huron	Erie	Ontario	Simcoe
Lake Whitefish
Superior		0.256 (35)	0.408 (35)	0.102 (31)	0.102 (35)	0.243 (35)
Michigan	0.23		<0.001 (37)	0.031 (31)	<0.001 (37)	0.031 (37)
Huron	0.16	0.77		0.113 (31)	0.001 (52)	0.021 (50)
Erie	−0.36	0.45	0.33		0.416 (31)	0.774 (31)
Ontario	0.33	0.66	0.51	0.16		0.350 (55)
Simcoe	0.24	0.41	0.39	−0.05	0.15
Cisco
Superior		0.903 (12)	0.715 (32)		0.452 (17)	0.903 (32)
Michigan	0.10		0.715 (12)		0.903 (11)	0.903 (12)
Huron	0.19	−0.35			0.715 (36)	0.715 (45)
Ontario	0.49	−0.06	0.20			0.858 (30)
Simcoe	0.02	−0.15	−0.14		−0.12

Statistically significant Pearson correlations are bold. The number of years used to calculate correlations are depicted in parentheses (above diagonal).

Lake Whitefish DFA model selection identified two models with ΔAICc ≤ 2, both of which included three latent trends and differed only in their error structure (Table 2). Both models produced near identical results, so we only interpreted results from the higher-ranked model. This model identified three latent trends across lakes using an error structure of the same variance with no covariance (Figs. 4 and 5). The first trend described a decline in YCS between 1960 and 1970, followed by a rapid sustained increase in YCS from 1970 to 1995, after which YCS steadily declined. All six lakes exhibited positive loadings with this trend, with the time series for lakes Michigan (factor loading = 0.31) and Huron (loading = 0.29) most strongly associated with this trend. The second trend showed high YCS between 1960 and 1970, followed by rapid declines reaching lows in the late 1970s, increases after 1980 to near average levels, a steady decline through the early 2000s, and a peak around 2010. Lakes Superior and Simcoe (loadings = 0.57 and 0.37, respectively) were strongly positively associated with this trend, whereas Lake Erie was strongly negatively associated (loading = −0.57). The third trend depicted lower than average YCS during the 1960s and 1970s followed by an increase in the 1980s–90s, a sharp decline around 2000, and a gradual increase thereafter. Lake Ontario was strongly positively associated with this trend (loading = 0.37). Lakes Michigan (loading = 0.11), Superior (loading = 0.07), and Huron (loading = 0.06) were also positively associated, while Lake Simcoe was negatively associated (loading = −0.11).

Figure 4.

Estimated latent trends (left panels) and corresponding factor loadings (right panels) for Lake Whitefish and Cisco year-class strength across study lakes based on dynamic factor analysis. Analyses were conducted separately for each species. Model selection supported three common trends for Lake Whitefish and one common trend for Cisco, respectively, using an error structure with the same variance with no covariance (Table 2). Factor loadings < 0.05 were not plotted. Y-axis limits are unitless and differ among rows.

Figure 5.

Fitted temporal trends (thick lines) in year-class strength (YCS) for each study lake (columns) based on the best supported dynamic factor model for Lake Whitefish (top row) and Cisco (bottom row), respectively. Shaded areas depict 95% confidence intervals. YCS estimates were z-score normalized prior to modeling (thin lines and points). Y-axis units are in standard deviations, where 0 (dashed line) denotes mean YCS relative to each model.

Table 2.

Akaike Information Criteria corrected for small sample sizes (AICc) model selection to identify the best supported error structure (variance-covariance) and number of underlying latent trends for dynamic factor analysis.

Variance-covariance structure	No. of trends	Log-likelihood	AICc	ΔAICc	AICc weight
Lake Whitefish
Same variance with no covariance	3	−293.63	621.42	0.00	0.67
Same variance and covariance	3	−293.44	623.31	1.90	0.26
Same variance with no covariance	4	−292.95	626.94	5.52	0.04
Same variance and covariance	4	−292.86	629.10	7.68	0.01
Same variance with no covariance	5	−292.10	629.92	8.51	0.01
Differing variance with no covariance	2	−298.00	632.43	11.02	0.00
Same variance and covariance	5	−293.08	634.26	12.84	0.00
Differing variance and covariance	3	−276.52	636.48	15.06	0.00
Same variance and covariance	2	−305.32	638.06	16.65	0.00
Same variance with no covariance	2	−306.43	638.08	16.66	0.00
Differing variance and covariance	4	−273.95	639.48	18.06	0.00
Differing variance and covariance	5	−273.28	643.66	22.24	0.00
Differing variance and covariance	2	−287.65	648.21	26.79	0.00
Differing variance with no covariance	1	−314.31	653.82	32.41	0.00
Differing variance and covariance	1	−299.54	659.34	37.92	0.00
Same variance and covariance	1	−327.72	671.98	50.57	0.00
Same variance with no covariance	1	−329.95	674.32	52.90	0.00
Cisco
Same variance with no covariance	1	-253.32	519.11	0.00	0.68
Differing variance and covariance	2	−233.16	521.87	2.76	0.17
Same variance and covariance	3	−246.61	523.70	4.59	0.07
Same variance with no covariance	3	−247.97	524.08	4.97	0.06
Same variance with no covariance	4	−247.97	528.79	9.68	0.01
Same variance with no covariance	2	−253.98	529.24	10.13	0.00
Same variance and covariance	2	−252.92	529.38	10.27	0.00
Differing variance and covariance	1	−242.77	530.69	11.58	0.00
Same variance and covariance	1	−258.79	532.22	13.11	0.00

Models that failed to converge after 1 000 000 iterations are not shown.

For Cisco, DFA model selection supported a single latent trend across lakes using an error structure of the same variance with no covariance (Table 2, Figs. 4 and 5). The trend was characterized by a peak in YCS in the early 1970s followed by a sharp decline to the late 1990s and an increase thereafter. Factor loadings with this trend were mostly positive and were highest for Lake Michigan (loading = 0.42), followed by Lake Simcoe (loading = 0.25) and Lake Huron (loading = 0.06), although Lake Michigan data coverage was limited to recent years (12 cohorts, 2004–2015). Lakes Superior and Ontario each had loadings very close to 0 (loadings = 0.01 and -0.01, respectively); consequently, Cisco YCS in lakes Superior and Ontario did not exhibit a trend that deviated from the long-term average through time (Fig. 5).

Discussion

We successfully estimated YCS for Lake Whitefish and Cisco cohorts spanning 1956–2015 across the Great Lakes and Lake Simcoe (Figs. 1 and 2). Our combined time series, comprising 38 long-term surveys with Lake Whitefish and Cisco catch and age data (Table A1), spans 1960–2019 and represents the most comprehensive dataset for analyzing Great Lakes coregonine recruitment compiled to date. We found that Lake Whitefish YCS was on average positively synchronous among lakes (Table 1), and Lake Whitefish in all six lakes experienced a period of elevated YCS during the 1980s–1990s (Figs. 4 and 5). In contrast, Cisco YCS was sporadic, highly variable around the long-term mean for each lake (Fig. 4 and 5), and not synchronous among lakes (Table 1). Our analyses also revealed that Lake Whitefish and Cisco YCS were not synchronous within lakes (Fig. 3). Collectively, our results suggest that the suite of factors regulating recruitment differs between these two species.

Our analyses revealed that Lake Whitefish exhibit positive synchrony in YCS across the Great Lakes and Lake Simcoe, and that YCS has declined in most of these lakes over the past two decades. Synchrony was most pronounced during the 1980s and 1990s when Lake Whitefish YCS was generally positive across all six lakes, coincident with the “phenomenal” recovery of Lake Whitefish stocks in the Great Lakes (Ebener 1997). This basin-wide recovery has been attributed to increased survival of multiple life stages resulting from favorable climatic conditions for embryos and larvae, declines in non-native pelagic planktivore populations (e.g. Alewife Alosa pseudoharengus and Rainbow Smelt Osmerus mordax), Sea Lamprey control, commercial harvest restrictions, and intermediate levels of primary productivity (Ebener 1997, Ludsin et al. 2001, Rook et al. 2022b). Although fishing pressure also declined during this period in Lake Simcoe (Dunlop et al. 2019), productivity remained relatively high (Evans et al. 1996, North et al. 2013), Rainbow Smelt populations increased (Evans and Waring 1987), and Sea Lamprey are absent (Docker et al. 2021). We therefore hypothesize that the observed spatial synchrony among Lake Simcoe and the Great Lakes reflects the importance of climatic processes operating across large spatial scales (Moran 1953, Koenig 2002) in regulating survival of Lake Whitefish early life stages (Brown et al. 2024). Specifically, cold winters are hypothesized to promote successful embryonic development (Brooke 1975), protect embryos from physical disturbance via ice cover (Freeberg et al. 1990), and depress abundances of predatory planktivores during larval emergence (Casselman et al. 1996). Lake Whitefish survival may be highest when cold winters are followed by warm springs (Taylor et al. 1987), which is thought to align larval emergence and the transition to exogenous feeding with the spring plankton bloom (match-mismatch hypothesis; Cushing 1990). Lake-specific factors are likely also important, as evidenced by our result that not all pairwise correlations among lakes were significant or positive (e.g. Lake Superior). Lake Whitefish YCS has declined in many of the Great Lakes in recent decades, although lake-specific trajectories varied. These declines were concurrent with reduced phosphorus loading, dreissenid mussel (Dreissena spp.) establishment, and loss of Diporeia, all of which are hypothesized to have acted in concert to limit prey availability, constrain carrying capacity, and exert strong density-dependent effects on Lake Whitefish (Rennie 2014, Gobin et al. 2015, Rook et al. 2022b). The Lake Superior ecosystem has largely been spared of these disruptions and Lake Whitefish populations remain stable (Fera et al. 2015, Carl et al. 2024); however, dreissenid mussel establishment (Trebitz et al. 2019) and warming winters (Mason et al. 2016) remain potential threats. Ultimately, a complex suite of interacting, context-dependent, and nonstationary factors may be driving spatial and temporal synchrony of Lake Whitefish YCS (Bourinet et al. 2023).

Cisco YCS was highly variable around the long-term mean in each lake and not synchronous among lakes. Dynamic factor analysis identified recent increasing trends in Cisco YCS in lakes Simcoe, Michigan, and, to a lesser extent, Huron. Populations in both Simcoe and Michigan are currently undergoing recovery from collapses that occurred by the 1980s–1990s (Finigan et al. 2018, Dunlop et al. 2019) and 1950s (Wells and McLain 1972, Claramunt et al. 2019), respectively. However, Cisco population abundances remain relatively low in both lakes and continued strong year-classes may be needed to rebuild spawning stocks. In contrast, Cisco YCS in lakes Superior and Ontario tracked the long-term mean through time, with no increasing or decreasing trends during our analysis frame. Large interannual fluctuations in YCS were observed across all lakes, providing additional evidence that sporadic recruitment may be intrinsic to Great Lakes Cisco populations (Stockwell et al. 2009). Our finding that Cisco YCS has been asynchronous among lakes was surprising, given that previous studies have found evidence of Cisco synchrony (Rook et al. 2012, Myers et al. 2015, Weidel et al. 2021) and that Lake Whitefish YCS was synchronous in this study. It is possible that low statistical power prevented detection of synchrony in our study, as Cisco YCS time series were generally shorter than those of Lake Whitefish due to data limitations. However, our study included multiple pairwise comparisons of 30+ year time series that still lacked evidence of synchrony. Rather, we hypothesize that our observed lack of synchrony could be due to lake-specific drivers (e.g. differences in food web structure; Gorman 2019, Brown et al. 2024) that outweigh regional climate effects (Myers et al. 1997). This could be why Myers et al. (2015) found that Cisco in Grand Traverse Bay (Lake Michigan) exhibited unique recruitment dynamics compared to other Upper Great Lakes populations. Populations could also be responding differently to similar climatic conditions, as Cisco are a notoriously plastic species (Eshenroder et al. 2016, George 2019) with evidence for differential spawning habitat preferences among lakes (Paufve et al. 2022). Embryonic incubation habitat varies along depth gradients among lakes, ranging from shallow shoals within protected embayments to deeper, offshore waters of Lake Superior. Consequently, embryos deposited at different depths might respond differently to similar climatic conditions during incubation, leading to asynchrony in recruitment. Future analyses of putative drivers of Cisco recruitment variability across the Great Lakes region could help test this hypothesis.

Our analyses highlight that Lake Whitefish and Cisco display fundamentally different recruitment dynamics in the Great Lakes region, despite sharing similar early life histories. The lack of synchrony between Lake Whitefish and Cisco YCS in our study aligns with previous findings that historical fishery-dependent CPUE for Lake Whitefish and Cisco in the Upper Great Lakes were largely independent of each other (Rook et al. 2021). These species-specific dynamics suggest that important recruitment drivers differ between species. Lake Whitefish and Cisco life histories are similar during early life stages (i.e. egg and larval; Brown et al. 2023) but diverge by the end of the larval stage when juvenile Lake Whitefish become demersal and benthivorous (Reckahn 1970) while Cisco remain pelagic (George 2019). This ontogenetic niche differentiation offers an opportunity to identify the life stages at which recruitment bottlenecks occur between species. Interspecific synchrony can be driven by climatic forcing (Liebhold et al. 2004), which would be expected to similarly affect early life stages of both species (Brown et al. 2024). Instead, divergent recruitment dynamics might arise during early life due to dissimilar responses to shared biophysical conditions, or during later life stages when Lake Whitefish and Cisco occupy different habitats with dissimilar biophysical conditions. Importantly, our approach to estimating YCS uses cohort-specific relative abundance starting at age-3 and integrates drivers of abundance from the time of egg deposition to vulnerability to sampling gears. Consequently, we hypothesize that the differing recruitment dynamics between species reflect differential drivers acting on post-larval life stages. Future research investigating recruitment drivers could deepen our understanding of the factors regulating these species across life stages (Brown et al. 2024).

The differing frequency of strong year-classes between Lake Whitefish and Cisco provides further evidence of species-specific recruitment dynamics in the Great Lakes region. Specifically, Cisco YCS was highly sporadic compared to that of Lake Whitefish, and this pattern was consistent among lakes for each species. The “boom-and-bust” recruitment dynamics observed for Cisco can result in populations being dominated by only a few strong year-classes (Rudstam et al. 1993, Yule et al. 2008), which can lead to uncertainty in fishery yields (Fisch et al. 2019) and may serve as an impediment to the reestablishment of extirpated populations (Rook et al. 2022a). In contrast, Lake Whitefish exhibited longer periods of sustained high or low recruitment. The universally above-average Lake Whitefish YCS observed during the 1980–90s fueled the growth of large-scale, economically valuable fisheries. However, this period of high YCS was relatively short-lived and it remains uncertain whether fishery yields during this period are a realistic expectation of potential stock productivity today or into the future, especially if unfavorable environmental conditions are responsible. Future research evaluating the role of intrinsic (e.g. life history traits; Winemiller 2005, Marquez et al. 2019) versus extrinsic (e.g. climatic forcing) factors in regulating each species’ population dynamics could help inform species- and lake-specific expectations for stock productivity.

In this study, we considered each lake a single population, but in many cases these data represent metapopulations of multiple spatially discrete subpopulations that may exhibit stock-specific recruitment dynamics (Ebener et al. 2021). For example, recent genetic evidence suggests that Lake Whitefish in Lake Michigan are spatially structured (Shi et al. 2022), and these populations have experienced differential recruitment trajectories in recent years (Modeling Subcommittee 2022). Stock-specific recruitment dynamics could be driven by local heterogeneity in biophysical conditions that relate to key drivers of recruitment (e.g. productivity and prey availability; Brown et al. 2024). Similarly, while we did not find evidence of interspecific synchrony at a lake-wide scale, it is possible that co-occurring spawning stocks of Lake Whitefish and Cisco could be more highly correlated at localized scales. However, although advanced genetic techniques (e.g. Euclide et al. 2022, Stott et al. 2022) and animal tracking technologies (e.g. Ebener et al. 2021, Gatch et al. 2023) have rapidly progressed understanding of contemporary stock structure for each species, spawning stock delineation remains a challenge in many areas of the Great Lakes. Though our lake-wide YCS estimates may have homogenized variation in trajectories of spawning populations within lakes, our approach allowed us to leverage the best available data to holistically describe YCS of Lake Whitefish and Cisco over multiple decades using a unified methodological framework. Future research aimed at quantifying the degree to which stock-specific recruitment dynamics vary across multiple spatial scales could increase understanding of metapopulation dynamics (e.g. portfolio effects; Schindler et al. 2010) and help optimize management decision-making for each species at relevant spatial scales (Ebener et al. 2021).

Our cross-lake synchrony results for each species disagreed with previous studies, which we suspect was due to differences in both the spatial scale of analysis and how recruitment was indexed among studies. While Lake Whitefish YCS was spatially synchronous in our study, Zischke et al. (2017) found that Lake Whitefish YCS—indexed using age-based catch-curves—was not synchronous among management units within lakes Superior, Michigan, and Huron. Importantly, the spatial scale of analysis should be taken into account when comparing inferences among studies of synchrony (Liebhold et al. 2004). It is possible that Lake Whitefish exhibit recruitment asynchrony among subpopulations within lakes (i.e. portfolio effects) and also exhibit broad-scale synchrony in response to concurrent, widespread ecosystem changes across lakes (e.g. dreissenid mussel establishment; Ebener et al. 2021). Previous studies of Cisco have documented recruitment synchrony at spatial scales spanning local to lake-wide (Rook et al. 2012, Myers et al. 2015, Weidel et al. 2021), but used methods to index recruitment (e.g. CPUE, catch-curve regression, stock-recruit models) that differed from our approach. These methodological contrasts are important to consider because how recruitment is characterized can impact resulting inferences, especially when the ages over which recruitment is indexed vary (e.g. Warren et al. 2024). For example, climatic forcing can act to synchronize abundances across populations during early life stages (e.g. eggs, larvae) but lake-specific processes might control survival at later life stages and ultimately desynchronize population dynamics (Dembkowski et al. 2016, Feiner et al. 2019). It is therefore possible that we may have detected stronger cross-lake recruitment synchrony for Cisco if our analyses focused on early life stages; however, abundances at those early life stages may not provide an accurate index of recruitment to adulthood if important recruitment bottlenecks occur during the juvenile life stage (Stige et al. 2013). To our knowledge, our study is the first to use the same methodological approach to compare recruitment dynamics between Lake Whitefish and Cisco or for either species across the entire Great Lakes basin. However, our ability to directly compare our results with previous studies is limited by the fact that this approach has not previously been used to investigate recruitment of either species in any of our study lakes. Future research could use our approach for YCS estimation at varying spatial scales to better understand how methodological decisions influence downstream inferences of recruitment synchrony.

We leveraged a diversity of long-term surveys to comprehensively describe Lake Whitefish and Cisco YCS in the Great Lakes region. While we attempted to account for differences in sampling design among surveys where possible, multiple sources of measurement error likely contributed to error in our estimates. Sampling gears that select for large-bodied fish (e.g. large-mesh gill nets) are typical for commercial monitoring surveys and surveys targeting Lake Trout (Salvelinus namaycush) but systematically observe older age classes of fish. While we were able to account for this by adjusting age ranges used to index YCS for each survey, the relative abundance of older age classes (i.e. those vulnerable to the fishery) is also influenced by fishing mortality (He et al. 2023). Importantly, our approach to YCS estimation assumes that cohort-specific relative abundance of likely recruits—age classes following the period of high natural mortality in early life—is a useful index for absolute recruitment to adulthood. Future research could evaluate this assumption by comparing our YCS estimates to time series of relative abundance of younger age classes; for example, existing surveys index abundance of age-0 Lake Whitefish in Lake Erie (Amidon et al. 2021) and age-1 Cisco in Lake Superior (Stockwell et al. 2009). Unfortunately, long-term monitoring programs that effectively index abundance of younger age classes of either species are rare across the Great Lakes, and efforts are ongoing to develop standardized sampling approaches (R. Tingley, U.S. Geological Survey, written comm., 18 Aug. 2023). We emphasize that this study leveraged a massive amount of data that has been collected by many different agencies over multiple decades. An analysis of this scale would not have been possible without these long-term survey programs and subsequent data management. Even so, availability of age data was our primary obstacle; we were unable to include many candidate surveys that reliably observed each species due to a lack of age data. Reliable age data is needed for demographic models that inform management and restoration decision-making (e.g. stock assessments, management strategy evaluations, structured population viability analyses), and continued cross-agency investment in age assessment across the Great Lakes basin may be needed to meet management and restoration goals. Lastly, we acknowledge that we did not propagate the uncertainty around our YCS estimates in the synchrony analyses; future research could build upon our study by developing statistical models of synchrony with error propagation and comparing the resulting inferences.

While there is a growing body of evidence that climatic processes operating across broad spatial scales act to synchronize recruitment dynamics of Great Lakes fishes, other biophysical processes that have been shown to regulate recruitment synchrony warrant further investigation. Climatic indices have been correlated with spatial synchrony in recruitment for Bloater (C. hoyi; Bunnell et al. 2010), Alewife (Warren et al. 2024), and Yellow Perch (Perca flavescens; Honsey et al. 2016) populations among the Great Lakes. Similarly, broad-scale climatic patterns are a well-documented driver of synchrony in marine fish population dynamics (Myers et al. 1997, Lehodey et al. 2006, Beaugrand et al. 2015). However, other physical and biological factors can be regionally correlated and interact with climate to influence synchrony (e.g. primary productivity and food web structure; Sheppard et al. 2019, Bourinet et al. 2023). Notably, Frank et al. (2016) found that fishing mortality could be an important driver of widespread spatial synchrony in Atlantic Cod (Gadus morhua) and cautioned that the prevailing view of climate as the predominant driver of synchrony resulted in other key drivers being overlooked. Alternatively, local conditions may supersede the synchronizing effects of climate (Marquez et al. 2023). For example, Vendace (Coregonus albula) recruitment was found to be spatially synchronous among lakes in Finland (Marjomäki et al. 2004) but not in Sweden (Axenrot and Degerman 2016), which Axenrot and Degerman (2016) hypothesized was due to interactions between lake-specific conditions (e.g. intraspecific competition among cohorts, fishing mortality) and regional climatic processes. Understanding of the processes regulating recruitment synchrony in the Great Lakes region could benefit from future studies considering a suite of climatic, physical, and biological variables in their analyses.

This study advances understanding of Lake Whitefish and Cisco recruitment dynamics in the Great Lakes region and highlights implications for their conservation and management. The spatially and temporally extensive time series of YCS reconstructed in this study enabled robust comparisons in recruitment variability among lakes and between species. Understanding the differences in the extent to which their recruitment is spatially synchronized helps optimize decision-making for each species at relevant spatial scales. For example, the lack of spatial synchrony for Cisco suggests that management and stewardship interventions may be most effective when tailored to the population characteristics and biophysical conditions of each lake, which is a guiding principle of ongoing Cisco restoration initiatives (Bunnell et al. 2023). Importantly, our YCS indices can be used to inform future research aimed at understanding the drivers of Lake Whitefish and Cisco recruitment variability at scales spanning local to basin-wide. Collectively, our results illustrate that the suite of factors regulating recruitment likely differ between Lake Whitefish and Cisco. Quantifying the important drivers of recruitment for each species, including the extent to which their relative importance varies among lakes, will be a key next step towards identifying which recruitment bottlenecks may be ameliorated by managers and stewards (e.g. biological control of invasive species, spawning habitat restoration) versus those that may serve as impediments to future recruitment (e.g. declining ice cover). Increased understanding of the mechanistic drivers of coregonine recruitment will be particularly beneficial in the face of accelerating ecosystem change across the range of these socioecologically valuable species (Austin and Colman 2007) and can help to mitigate uncertainty regarding how multiple populations will respond to future conditions (Mulvaney et al. 2014, DeVanna Fussell et al. 2016).

Acknowledgments

We are indebted to the countless biologists, vessel crews, and support staff who collected, processed, and maintained decades of survey data. We specifically thank the following agencies and individuals who shared data for use in this project: Bay Mills Indian Community (Jack Tuomikoski); Grand Traverse Band of Ottawa and Chippewa Indians; Great Lakes Fishery Commission Lamprey Database and representatives of the Lake Michigan Technical Committee (Dave Clapp, Ben Dickinson, Kevin Donner, Chuck Madenjian, Dan Makauskas, Archie Martell, Cheryl Masterson, Rebecca Redman, Ben Turschak, Susan Wells); Little Traverse Bay Bands of Odawa Indians (Kevin Donner, Gary Michaud, Nicholas Moorman, Katherine Skubik); Michigan Department of Natural Resources (Tracy Claramunt, David Fielder); Minnesota Department of Natural Resources (Cory Goldsworthy); Ontario Ministry of Natural Resources and Forestry (Megan Belore, Chris Davis, Rachel Henderson, Blair Wasylenko, Zoe Zrini); Red Cliff Band of Lake Superior Chippewa (Ian Harding); Sault Ste Marie Tribe of Chippewa Indians (Brad Silet, Amanda Stoneman); U.S. Fish and Wildlife Service (Jose Bonilla Gomez, Kevin McDonnell); and Wisconsin Department of Natural Resources. We also thank Nicole Berry, two anonymous reviewers, and journal staff who helped improve this manuscript. This is NOAA GLERL contribution number 2055. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Author contributions

All authors contributed to the development of this manuscript. TAB, LGR, SAS, DBB, BCW, and AEH conceptualized the study. Funding was acquired by TAB, LGR, SAS, PR, JBS, TJT, EO, JXH, JLJ, EKB, EB, SJ, SAP, JT, JAS, SPH, DBB, BCW, and AEH. Data collection and curation was led by TAB, PR, JBS, TJT, CH, EO, JXH, JLJ, JEB, SJHB, EKB, EB, AC, SJ, JT, DDC, and SPH. TAB led the study design and methodology, formal analysis, and writing of the original draft with supervision from LGR, SAS, DBB, BCW, and AEH and guidance from all authors. All authors critically reviewed and edited the manuscript.

Funding

This work was supported by the Great Lakes Fishery Commission (2022_HON_441017).

Data availability

The derived YCS indices generated in this study will be shared on reasonable request to the corresponding author. The long-term survey data underlying the YCS indices were provided by each respective agency to the corresponding author and cannot be shared publicly per data sharing agreements. Data can be obtained by request to the agency responsible for each survey (see Supplementary Material S1).

Appendix

Table A1.

Surveys used to estimate year-class strength by lake, species, and sampling gears.

		Lake Whitefish		Cisco
Survey Name	Type	Sampling gears	Year-classes included in analysis	Sampling gears	Year-classes included in analysis
Lake Superior
BMIC Commercial Fisheries Monitoring	FDS	GN-LG; TN	1986–2014	−	−
BMIC Lake Trout Assessment	FIS	GN-LG	1993–1999; 2000–2013	−	−
BMIC Lake Whitefish Assessment	FIS	GN-GRD	2000–2015	GN-GRD	2000–2005; 2008–2015
MDNR Statewide Angler Survey Program	FDS	ANG	1981–2011	−	−
MNDNR Lake Superior Cisco Assessment	FIS	−	−	GN-SM	1994–1996; 2005–2015
MNDNR Lake Superior Spring and Fall Cisco Commercial Sampling	FDS	−	−	GN-SM	2005–2015
MNRF Lake Superior Commercial Catch Sampling	FDS	GN-LG	1990–2012	GN-SM	1991–2015
MNRF Lake Superior Fish Community Index Netting	FIS	GN-GRD	2004–2015	GN-GRD	2004–2015
MNRF Spawning Cisco Hydroacoustic Survey	FIS	−	−	GN-GRD	2008–2015
RC Fall Spawning Lake Whitefish Survey	FIS	GN-LG	1981–2010	−	−
ST Commercial Harvest Monitoring	FDS	GN-LG; TN	2003–2013	−	−
WDNR Lake Superior Commercial Fishery Monitoring	FDS	GN-LG	1997–2013	GN-SM	1997–2015
WDNR Lake Superior Fall Spawning Cisco Survey	FIS	−	−	GN-GRD	1984–2006; 2011–2015
WDNR Lake Superior Spring Lake Trout Assessment	FIS	GN-LG	1981–1986; 1992–2013	−	−
WDNR Lake Superior WI-1 and WI-2 Summer Graded Mesh Surveys	FIS	GN-GRD	1993–2015	GN-GRD	1997–2015
Lake Michigan
GTB Commercial Harvest Monitoring	FDS	GN-LG; TN	2008–2014	−	−
LTBB Commercial Harvest Monitoring	FDS	GN-LG; TN	1993–2013	−	−
MDNR Statewide Angler Survey Program	FDS	ANG	1979–2014	ANG	2005–2015
MDNR Northern Lake Michigan Reef Gillnetting Surveys	FIS	GN-GRD	2004–2009; 2011–2015	GN-GRD	2004–2014
LMTC Lakewide Assessment Plan and Fishery-Independent Whitefish Survey	FIS	GN-GRD	1993–2015	GN-GRD	2004–2015
ST Commercial Harvest Monitoring	FDS	GN-LG; TN	2003–2014	−	−
ST Fishery-Independent Whitefish Population Assessment	FIS	GN-LG	1994–2013	−	−
WDNR Lake Michigan Dockside Commercial Whitefish Fishery Monitors	FDS	GN-LG; TN	1983–2014	−	−
WDNR Lake Michigan Juvenile Lake Whitefish Population Assessments	FIS	GN-GRD	1985–2015	−	−
WDNR Lake Michigan Adult Lake Whitefish Population Assessments	FIS	GN-GRD; EF	1985–2015	−	−
Lake Huron
MDNR Statewide Angler Survey Program	FDS	−	−	ANG	1994–2011
MDNR Lake Huron Commercial Fishery Biosampling	FDS	TN	1975–2015	−	−
MDNR Lake Huron Annual Spring Gillnetting Surveys	FIS	GN-GRD	1965–1984; 1991–2015	GN-GRD	1969–1980; 1995–1997; 2010–2015
MDNR Les Cheneaux Islands Fish Community Survey	FIS	−	−	GN-GRD	1965–1974; 1996–2002; 2006–2011
MNRF Lake Huron Commercial Catch Sampling	FDS	GN-LG; GN-SM; TN	1969–2015	GN-LG; GN-SM; TN	1968–1982; 1986–2009
MNRF Lake Huron Offshore Index Assessment	FIS	GN-GRD	1964–2015	GN-GRD	1973–2015
ST Drummond Island Refuge Assessment–Lake Herring	FIS	−	−	GN-LG; GN-SM	2004–2014
ST Commercial Harvest Monitoring	FDS	GN-LG; TN	2004–2015	GN-LG; TN	2003–2013
ST Drummond Island Refuge Assessment–Lake Trout	FIS	GN-SM	1993–2015	GN-SM	1991–2015
ST Fishery-Independent Whitefish Population Assessment	FIS	GN-LG	1995–2015	GN-LG	2003–2013
USFWS Lake Huron Lake Whitefish Adult Assessment Survey	FIS	GN-GRD	1997–2015	−	−
Lake Erie
MNRF Lake Erie Commercial Catch Sampling	FDS	GN-LG	1990–2015	−	−
MNRF Lake Erie Partnership Survey	FIS	GN-GRD	1985–2015	−	−
Lake Ontario
MNRF Lake Ontario and Bay of Quinte Fish Community Index Gillnetting	FIS	GN-GRD	1956–2015	GN-GRD	1957–1984; 1988–1990; 2005–2015
MNRF Lake Ontario Commercial Catch Sampling	FDS	GN-LG; TN	1987–2013	TN	2000–2001; 2005–2015
Lake Simcoe
MNRF Fall Index Trap Netting	FIS	TN	1958–1962; 1966–2000; 2002–2004	TN	1971–2000
MNRF Offshore Benthic Index Netting	FIS	GN-GRD	1998–2015	GN-GRD	1998–2015
MNRF Lake Simcoe Winter Creel Survey and Winter Angler Catch Sampling	FDS	ANG	1959–1961; 1966; 1969–1977; 1980; 1982–1994; 1996–2004; 2006–2008	ANG	1971–1995

Full descriptions of each survey are available in Supplementary Material S1. Lakes include Superior (SU), Michigan (MI), Huron (HU), Erie (ER), Ontario (ER), and Simcoe (SC). Surveys are categorized as either fishery-independent (FIS) or fishery-dependent (FDS) survey types. Sampling gears include graded-mesh gill nets (GN-GRD), small-mesh gill nets (GN-SM), large-mesh gill nets (GN-LG), trap nets and other impoundment nets (TN), electrofishing (EF), and angling (ANG). Agency names are acronymized: Bay Mills Indian Community (BMIC); Grand Traverse Band of Ottawa and Chippewa Indians (GTB); Lake Michigan Technical Committee (LMTC); Little Traverse Bay Bands of Odawa Indians (LTBB); Michigan Department of Natural Resources (MDNR); Minnesota Department of Natural Resources (MNDNR); Ontario Ministry of Natural Resources and Forestry (MNRF); Red Cliff Band of Lake Superior Chippewa (RC); Sault Ste Marie Tribe of Chippewa Indians (ST); U.S. Fish and Wildlife Service (USFWS); Wisconsin Department of Natural Resources (WDNR).

Table A2.

Age classes used to index year-class strength by lake, sampling gear, and species.

		Age Ranges
Lake	Sampling Gear	Lake Whitefish	Cisco
Superior	ANG	3–5	−
	GN-GRD	3–5	3–5
	GN-LG	6–7	−
	GN-SM	−	4–6
	TN	5–6	−
Michigan	ANG	5–6	4–6
	EF	6–8	−
	GN-GRD	3–5	3–5
	GN-LG	6–7	−
	TN	5–6	−
Huron	ANG	−	3–5
	GN-GRD	3–5	3–5
	GN-LG	4–6	6–7
	GN-SM	3–4	4–6
	TN	4–6	6–8
Erie	GN-GRD	3–5	−
	GN-LG	4–6	−
Ontario	GN_TN	5–7	−
	GN-GRD	3–5	3–5
	GN-LG	6–8	−
	TN	6–8	3–5
Simcoe	ANG	6–8	5–7
	GN-GRD	3–5	3–5
	TN	6–8	5–7

Sampling gears are abbreviated: graded-mesh gill nets (GN-GRD); large-mesh gill nets (GN-LG); small-mesh gill nets (GN-SM); trap nets and other impoundment nets (TN); electrofishing (EF); angling (ANG).