Landscape and biophysical controls of lake productivity to inform evaluation of sockeye salmon (Oncorhynchus nerka) populations in data‐limited regions

Landscape models are increasingly used to classify and predict the structure and productivity of data‐limited aquatic ecosystems. One such suite of ecosystems is on the remote North and Central Coast (NCC) of British Columbia, where sockeye salmon (Oncorhynchus nerka) rear in more than 150 lakes. Given their remoteness and limited resources for assessment, limnological and population monitoring in many of these lakes has been periodic or absent, limiting understanding of the status of populations and their habitats. Lake photosynthetic rate (PR) estimates are foundational to models of sockeye salmon nursery lake productive capacity. Using data from 61 lakes across the NCC, we compared a suite of landscape and lake variables in an information theoretic framework producing a set of models relating these characteristics to lake PR. A categorical variable related to lake biogeochemistry—whether a lake is humic stained, clear, or glacially turbid—was the most important variable predicting lake PR and was included in all models. Lake surface area relative to upstream catchment size and lake perimeter‐to‐surface‐area ratio were also important, with smaller upstream catchments yielding higher production, and high shoreline complexity correlated with lower productivity as measured by limnetic PR. Model‐averaged predictions of PR from the four models with the lowest residual error were created for 96 lakes currently lacking limnological assessments. These landscape models represent a valuable starting point for evaluating lake‐specific carrying capacities for data‐poor sockeye salmon populations under Canada's Wild Salmon Policy.

A key applied challenge for ecologists is to quantify the productivity and capacity of ecosystems for the management and conservation of species in data-limited regions. In light of this challenge, landscape-scale models have been proposed as a tool to leverage information on regional and local habitat conditions for prediction and prioritization of conservation and management over broad spatial scales (Soranno et al. 2010;Schwenk and Donovan 2011). Ecosystems are governed by an interacting hierarchy of physical and biological processes, ranging from local habitat conditions and species interactions to broad regional differences in climate, geology, and community structure (Levin 1992;Fergus et al. 2011). Thus, proximal habitats often share similar geomorphic and habitat conditions and exhibit similar patterns of community productivity and structure (Legendre 1993;Lichstein et al. 2002). The effects of landscape and climate conditions can therefore be modeled across systems by drawing on highquality information from a representative sample of locations to build broader regional understanding of ecosystem conditions (Legendre and Fortin 1989;Turner and Gardner 2015).
In recent years, there has been a dramatic increase in the use of geospatial data paired with data on species distribution or abundance to evaluate habitat suitability and use. By understanding the underlying landscape drivers of species distribution or abundance, ecosystem productivity, and trophic structure (Elith and Leathwick 2009;Soranno et al. 2010), these models can support conservation planning and prioritization, and extrapolate this inference to unsampled or remote habitats. For example, spatial models of mountain caribou habitat that account for both land cover and geomorphology have shown promise for prioritizing conservation and land use planning (Johnson et al. 2004). Indeed, models of species distribution data in relation to landscape level variation in geomorphology and habitat structure are used to predict the distribution and abundance of a wide variety of at-risk species, ranging from fish to fishers (Carroll et al. 1999;Pess et al. 2002).
In lake ecosystems, productivity and community structure are driven by a complex suite of physical and chemical processes, and biological interactions such as predation and competition (Goldman and Horne 1983;Carpenter and Kitchell 1988). The hydrology and geomorphology of lake catchments interact with regional-scale drivers such as climate and land use to govern water chemistry and nutrient availability, flushing rates, temperature, and euphotic depth (Kratz et al. 1997;Kamenik et al. 2001;Fergus et al. 2011). Through these diverse pathways, landscapes contribute to the regulation of primary productivity and the abundance of species at higher trophic levels within lakes (Hershey et al. 1999;Quinlan et al. 2003).
Limnologists and fisheries biologists have long sought to understand linkages between lake ecosystem properties and fisheries yields (Northcote and Larkin 1956;Hanson and Leggett 1982;Jones and Hoyer 1982). Early efforts examined the relationship between fisheries productivity and simple metrics such as total dissolved solids and depth (Northcote and Larkin 1956;Ryder 1965). More recently researchers have measured the rate of primary production within the limnetic food web of lakes and found that photosynthetic rates (PRs) are highly correlated with fish biomass (McConnell et al. 1977;Downing et al. 1990;Hume et al. 1996). While predictive models of lake productivity, in isolation, necessarily simplify ecosystem dynamics, the similarities among lakes and regions provide the opportunity to make generalizable predictions about ecosystem conditions and processes across broad spatial scales (Wagner and Schliep 2018). Landscape models of ecosystem productivity are therefore attractive for managers seeking to inform management of fisheries with limited resources.
Across landscapes, watersheds integrate catchment-scale variation in hydrology and geomorphology, as well as regional forcings including climate, biogeography, and anthropogenic disturbance, producing regionally coherent patterns in trophic structure, organic material and nutrient delivery, and productivity (Fergus et al. 2011;King et al. 2019). On the Pacific coast of Alaska and British Columbia, lakes within the same biogeographic regions often share similar water chemistry, hydrology and climate, resulting in broad, regionally defined patterns of lake productivity (Stockner and MacIsaac 1996). For the purposes of management and prediction, lakes in coastal British Columbia and Alaska are commonly classified as either stained, clear, or glacially turbid (Edmundson and Mazumder 2001;Shortreed et al. 2007). Lakes in these groupings are spatially clustered, reflecting regional gradients in hydrologic and geomorphic conditions. Low elevation coastal lakes are typically stained from high concentrations of dissolved organic materials (DOM), clear water lakes are primarily found interior watersheds with lower precipitation and lower rates of DOM delivery, and glacial lakes are found in higher elevation catchments. Consequently, whether a lake is stained, glacial or clear is related to similarities in the physical, chemical, and biological conditions influencing primary, secondary, and fish production (Lloyd et al. 1987;Stockner and MacIsaac 1996;Shortreed et al. 2001). Therefore, previous studies have used these three categories of water clarity-clear, glacial, and stained-when seeking to characterize patterns of sockeye salmon nursery lake productivity (Edmundson and Carlson 1998;Shortreed et al. 2007).
Sockeye salmon (Oncorhynchus nerka) typically depend upon lake habitats for juvenile rearing. They have a pan-North Pacific range, and are of major cultural, economic, and ecological significance. Thus, models linking lake productivity to landscape and lake conditions could improve conservation and management prospects for many data-limited populations. Due to the tight coupling of juvenile sockeye salmon production with the pelagic food webs of their natal rearing lakes, their populations are often limited by lake productivity and size (Juday et al. 1932;Shortreed et al. 2001). In recent decades, researchers and managers in Alaska and British Columbia have developed limnological rearing capacity models for sockeye salmon nursery lakes (Koenings and Burkett 1987;Hume et al. 1996;Shortreed et al. 2000). In British Columbia, data on lake PRs have been used to predict juvenile rearing capacity for sockeye nursery lakes Shortreed et al. 2000). Model outputs have been used to inform stock conservation and harvest strategies, and to estimate stock-specific recovery potentials. For example, informative priors based on lake rearing capacity are often used in Bayesian stock assessment to improve model fits and reduce uncertainty in estimates of biological or management benchmarks (Grant et al. 2011). These approaches are particularly relevant in populations where stock-recruit data are scarce (Cox-Rogers et al. 2010).
Linking landscape conditions to lake productive capacity would be particularly useful in remote regions where population monitoring and full growing season PR data are challenging and costly to obtain. On the North and Central coast (NCC) of British Columbia, sockeye salmon support important subsistence and commercial fisheries, but many populations have shown declining productivity and abundance in recent decades (Peterman and Dorner 2012). There are at least 157 NCC lake systems supporting populations of sockeye, and each population is considered sufficiently genetically or demographically distinct to justify protection and management as a conservation unit (Holtby and Ciruna 2007). From low elevation NCC coastal lakes to the mountainous watersheds in the interior, such as those of the Skeena and Nass Rivers, these sockeye salmon populations represent the diverse evolutionary Atlas et al.
Models of sockeye lake rearing capacity and ecological legacy of their species (Wood et al. 1994). Given the remote nature of these lakes and the limited resources for fish and fish habitat monitoring, many stocks lack basic assessment information. To understand the physical, chemical, and biological factors limiting freshwater productivity of sockeye salmon in Canada, Fisheries and Oceans Canada's (DFO) Lakes Research Program conducts limnological surveys of sockeye salmon nursery lakes to estimate lake food web productivity and structure, producing habitat-based estimates of optimal adult abundance and juvenile production through application of the PR Model Shortreed et al. 2000;Cox-Rogers et al. 2010). For this purpose, growing-season (i.e., May to October) limnological assessments have occurred for 61 lakes across the NCC (Shortreed et al. 2000(Shortreed et al. , 2007. In the NCC, these physical, chemical, and biological conditions have not been integrated into landscape-scale predictive models to provide broader insights into regional patterns of PR and sockeye salmon production. Therefore, there is interest in developing a regional-scale predictive model which links easily derived landscape metrics to PR, providing predictions of productive capacity for the 96 lakes where data are currently lacking. Our goal was to predict PR across the NCC region by quantifying the linkages between biophysical and geomorphic conditions at the landscape level and lake productivity. Using models linking landscape and hydrological conditions in 61 lakes to their observed PR, we developed a quantitative framework for predicting PR across 96 previously unsampled NCC nursery lakes. These estimates can inform future detailed assessment, management planning, and evaluation of stock status, particularly when habitat-based productive capacity estimates are used as benchmarks for population assessment in data-limited sockeye populations (Cox-Rogers et al. 2010). More broadly, this work highlights opportunities to apply landscape approaches to inform the management and conservation of data-limited and culturally important ecosystems.

Limnological data and PR estimates
Estimates of individual lake PR were obtained from a series of DFO studies estimating annual primary productivity and juvenile rearing capacity for each of the 61 lakes (Supporting Information Table S1), falling within 49 watersheds. These data were collected over the span of 30 yr between 1978 and 2008. A few lakes (n = 7) had 2 yr of data and Meziadin was sampled three times bringing the total number of data points to 70. Lakes did not exhibit systematic patterns of directional changes in PR so values were averaged for these lakes for the purpose of analysis (Supporting Information Fig. S2).
PR estimates were derived using in situ light and dark bottle incubations, with samples spanning the euphotic zone of each lake, from which the autotrophic uptake of 14 C isotopes was measured (Shortreed et al. 2000). The number of sampling locations varied from one to five across the 61 study lakes, depending on lake size and mixing. At each location, the euphotic zone was sampled at 4-6 depths based on thermocline and compensation depths. Water from each depth was collected in three 125-mL clear and two 125-mL opaque bottles. Bottles were incubated at their respective depth for 2 h between 09:00 and 12:00 h, and then immediately transported for lab analysis in light proof boxes. A 40-mL aliquot of this water was then vacuum filtered on filters with 0.2 and 2.0 μm pore sizes, then dried before being placed into scintillation vials. Shortly thereafter samples were counted on scintillation counter  to estimate hourly PRs. These measurements were integrated with concurrent data on euphotic depth and morphometry to yield daily PR estimates (mg C m −2 d −1 ). In most instances, lakes were monitored monthly throughout the growing season (e.g., May-October), and growing season averages were used. However, in more remote lakes only a single late-summer sampling occurred, and seasonal mean PR estimates were estimated from the equation for NCC lakes (PR seasonal mean = 0.7479[PR fall ]; r 2 = 0.60, p < 0.05, n = 113) (Cox-Rogers et al. 2004).
In the sockeye salmon PR rearing capacity model developed by Hume et al. (1996), estimates of mean growing season production are converted to total lake-wide growing season production (PR TOTAL ) by multiplying mean daily PR by lake surface area (m 2 ) of the focal lake and growing season length (183 d). Total growing season production is related to sockeye salmon juvenile rearing capacity because juvenile production is often limited by the productivity of lake pelagic ecosystems Shortreed et al. 2000). This information is then converted into an estimate of sockeye spawner abundance at carrying capacity (S MAX_PR ), using a model presented in Hume et al. (1996) that multiples total annual photosynthetic production or PR TOTAL (tons C yr −1 ) by the constant 187 (spawners tons C −1 ) to yield an estimate of the number of adult spawners required to maximize smolt production (Eq. 1). Given the scarcity of robust time series of sockeye spawner abundance, these PR model estimates of lake capacity have been used as stock-specific benchmarks for evaluation of stock status in data-limited sockeye populations in the Skeena watershed (Cox-Rogers et al. 2010).
Landscape variables Previous investigations of sockeye salmon nursery lake productivities have revealed major biological and physical differences related to lake water clarity (i.e., stained, glacial, clear) (Edmundson and Carlson 1998;Edmundson and Mazumder 2001;Shortreed et al. 2007). In general, these lake clarity groupings are useful, because they are associated with important differences in hydrology, geomorphology, climate, and land cover. While lake clarity assignments for sockeye rearing lakes have traditionally been made qualitatively (Shortreed Atlas et al. Models of sockeye lake rearing capacity et al. 2007), these lake clarity groups exhibit different ratios of euphotic zone depth and Secchi depth, owing to differences in light scattering and absorption between clear, glacial, and stained lakes (Koenings and Edmundson 1991;Shortreed et al. 2007). Preliminary data exploration revealed strong gradients in lake productivity and biophysical conditions associated with these water clarity groups, and we included lake water clarity as a categorical, independent variable in all models. Since no limnological data were available for the 96 unsampled lakes, we assigned a water clarity based on the values assigned to surrounding lakes, where lakes clustered in a Principal Components Analysis biplot (Supporting Information Fig. S1), and then visually confirmed lake color using Google Earth™ imagery as well as expert knowledge from DFO limnology and sockeye stock assessment biologists (D. Selbie and S. Cox-Rogers).
We selected a suite of 11 landscape and lake variables, derived using ArcMap 10.3 Geographic Information Systems (GIS) software, to be used as potential predictor variables for lake productivity. These variables were chosen based on specific hypothesized relationships between landscape and productivity of freshwater ecosystems (Table 1) including: watershed area upstream from the outlet of each respective sockeye-bearing lake (n = 157), watershed elevation (mean and max), mean watershed slope, the proportion of upstream watershed occupied by lakes and bogs (upstream lake), the proportion of upstream watershed occupied by the focal lake (lake-to-watershed proportion), % glacial cover, % forested, annual precipitation, and mean growing season temperatures. For climatic variables, the growing season was defined as May through August. We obtained historical climate data using the ClimateBC tool (Wang et al. 2016), and estimated annual precipitation and mean growing season temperature across each watershed for the climate normal period (1975-present), overlapping the period when PR data were collected. Watershed area, elevation and slope were derived using a 20 m digital elevation model. Glacier area was obtained from the Randolf Table 1. Predicted relationship between 11 candidate variables and lake PR.

Landscape feature Hypothesis
Latitude Latitude influences temperature and growing season length, and is correlated with lake productivity (Håkonsen and Boulion 2001).
Water clarity category (clear, glacial, stained) Humic stained, glacially turbid, and clear water lakes exhibit distinct physical and biological conditions. These differences in nutrients, light, pH, seasonal temperature stratification, and trophic structure act to control primary productivity contributing to regional variation in lake productivity (Lloyd et al. 1987;Stockner and MacIsaac 1996;Jansson et al. 2000). Distance from coast Climatic and biogeographic variation from coastal to interior watersheds creates gradients in precipitation, temperature, land cover, hydrology, nutrient, and sediment delivery. Accordingly, productivity varies moving inland from low coastal watersheds to mountainous and interior plateau lakes (Shortreed et al. 2007).
Mean watershed elevation Elevation mediates temperature and growing season length, snowpack, and hydrology (Isaak and Hubert 2001;Lisi et al. 2015). These physical factors underpin rates of biological productivity including lake PR.

Maximum watershed elevation
Maximum elevation captures the degree to which snowpack contributes to discharge during the growing season, stabilizing temperature, and controlling the timing of water and nutrient delivery (Lisi et al. 2015).
Upstream lake (upstream lake area/ drainage area) The amount of upstream lake and wetland can influence water chemistry, temperature, nutrient availability, primary and secondary production (Quinlan et al. 2003;Fergus et al. 2011;Sadro et al. 2012).
Perimeter-to-area ratio Lake morphometry is related to primary production (Oglesby 1977). High lake perimeter ratio indicates greater extent of littoral habitat, increasing coupling between littoral and terrestrial habitat with lake nutrient dynamics and food webs (Vadeboncouer et al. 2002).
Watershed slope Watershed slope is related to peak discharge and flushing rates, control nutrient delivery, and export (Kamenik et al. 2001). Low watershed slopes may be associated with increased temperature accumulation (Lisi et al. 2015) and high delivery of dissolved organic carbon (Rasmussen et al. 1989).
Lake-to-watershed proportion (lake area/drainage area) Lakes with larger upstream watersheds receive higher contributions of water and organic matter from upstream catchment (Rasmussen et al. 1989). Large upstream drainages relative to lake area influence pathways of water delivery and flushing rates, nutrient delivery, and retention (Kratz et al. 1997).
Watershed % tree cover Tree cover is related to temperature (Isaak and Hubert 2001), weathering, and nutrient delivery (Kamenik et al. 2001) in lotic ecosystems, influencing rates of primary production. Watershed % glacier Glaciers linked to variation in timing and intensity of discharge, temperature, as well as sediment and nutrient delivery. Glacial turbidity may also limit euphotic depth and hinder primary productivity (Lloyd et al. 1987).

Atlas et al. Models of sockeye lake rearing capacity
Glacier Inventory (RGI Consortium 2017), vegetation data (tree cover specifically) was obtained from the Vegetation Resources Inventory with data inputs ranging from 1990 to 2018, and percent cover data were calculated using respective watershed area. Lake specific information including perimeterto-area ratio and distance from coast was obtained from the British Columbia Freshwater Atlas, with the distance to coast (m) variable being measured as a straight line from the lake outlet to the nearest coastline. Because of low topographical relief, we were not able to obtain separate watershed characteristics for the lower and middle Mikado lakes, and therefore combined the two lakes into a single data point for analysis. Four variables, the proportion of watershed area occupied by the sockeye-rearing lake (referred to hereafter as "lake-towatershed proportion"~lake area/watershed area), the proportion of upstream watershed area occupied by lakes and bogs (referred as "upstream lake"~upstream lake area/watershed area), and the landcover variables (proportion glacial cover and tree cover) were transformed using logit transformations commonly applied to proportion data. Another derived variable (termed "perimeter-to-area ratio") intended to capture the littoral contribution to lake surface area, was estimated as the ratio of perimeter (m) to surface area (m 2 ) for each lake. We naturallog transformed these ratio data.

Model selection and lake predictions
Preliminary data exploration revealed a high degree of correlation among some landscape variables. To reduce the number of candidate variables and eliminate problems associated with collinearity, we performed stepwise variable reduction by estimating generalized variance inflation factors (VIF) for each continuous landscape variable among the candidate set and sequentially eliminated those with the highest VIF scores greater than 10 (Craney and Surles 2002). This procedure wasrepeated until only variables with a VIF lower than 10 remained among the candidate set, eliminating mean watershed elevation, max watershed elevation, and distance from coast from our candidate variable set. We then compared a suite of linear mixed-effects models relating the remaining landscape and climate variables (1. latitude, 2. water clarity category, 3. upstream lake, 4. perimeter-to-area ratio, 5. watershed slope, 6. lake-towatershed proportion, 7. tree cover, 8. glacier cover, 9. growing season temperature, 10. annual precipitation) to lake PR for the 61 lakes with available PR data. All models included a random intercept (μ j ) based on their effect of watershed identity to account for the fact that multiple lakes were nested within some watersheds (j) (e.g., Atnarko, Babine, Banks, Kispiox, Mikado, Zymoetz, and other sub-basins in the Skeena watershed with multiple rearing lakes). Lake PR was natural-log transformed to meet the assumptions of normality associated with linear modeling. We limited the number of interactions considered to two potential interactions: one between our categorical variable water clarity and latitude, and another between water clarity and growing season temperature. Given differences in light penetration, heat retention, and nutrient availability, we hypothesized that changes in temperature and growing season length associated with latitude could manifest differently across water clarity categories. Furthermore, data visualization suggested differences in the slope of the relationship between latitude and productivity across the lake clarity types.
Given the high number of potential combinations of landscape and climate variables, we winnowed our candidate model set using Akaike's Information Criterion for small sample sizes (AICc) producing a set of models that best fit the data without overfitting (Burnham and Anderson 2002). We evaluated all additive combinations of variables, using the dredge function in the R package MuMIn (Barto n 2018). Because AICc identified several models with a high degree of support, we ranked models based on their δAICc score-the difference between the AICc score of a given model and the top ranked model. A total of 19 models fell within this threshold (Burnham and Anderson 2002).
To compare goodness of fit among these models, we computed root-mean-square-error (RMSE) for each model. We then generated predictions of PR for the 96 unsampled lakes using the four models with the lowest RMSE using the predict function in R (R-Core Development Team). Mixed-effects models pose unique mathematical challenges for prediction since random effects and fixed effects are estimated as separate variance components (Gelman and Hill 2007), and estimating error for random effects requires bootstrapping or Bayesian MCMC. Confidence intervals (CIs) for predicted values were therefore generated using standard errors from the model fit without random effects, reflecting uncertainty at the fixedeffect level. Only eight of the 96 unsampled lakes shared watersheds with those in our 61 lake sample. We therefore predicted PR for these 96 lakes using coefficient values and uncertainty from fixed effects, and incorporated random intercepts post hoc for lakes within sampled watersheds. To eliminate bias associated with backtransforming from the natural log scale, we performed a correction following the method recommended by Dambolena et al. (2009), computing the mean square error (ε) from the residual error in the fit of our model, and multiplying predicted PR by e 0.5ε .
To evaluate goodness of fit, we performed leave-one-out cross validation for the four models with the lowest RMSE to estimate out-of-sample predictive performance. Data points were sequentially excluded and the model refit, and used to predicted PR for the missing lake.
Mean and 95% confidence intervals of these PR predictions were compared with observed values for each lake in the sample (Supporting Information Fig. S3, Table S2).

Variation in lake conditions and PR
Geomorphology, hydrology, and climate vary widely among sockeye lakes in the NCC. Lake elevations ranged Atlas et al.
Models of sockeye lake rearing capacity from a minimum of 5 m above sea level (ASL) in Bonilla, Curtis, and Moore Lakes, to 1448 m ASL for Johanson Lake in the Sustut watershed, a tributary of the Skeena River. Accordingly, annual rainfall and temperature varied dramatically between low elevation coastal and mountainous interior watersheds. Mean annual precipitation in coastal watersheds was 3879 mm, compared to 974 mm in interior watersheds. Growing season temperatures also reflected strong climatic gradients, with interior watersheds experiencing a mean air temperature of 9.21 C from May through August, and coastal watersheds having a mean air temperature of 11.65 C during that same period. However, lake catchment and climate characteristics were roughly equivalent among sampled and unsampled lakes in our study, allowing predictions of PR in unsampled lakes from coefficients estimated in systems with PR data (Supporting Information Table S3). Lake clarity category showed strong spatial structure and was associated with a suite of limnological characteristics. Stained lakes tended to be concentrated along the coast, while clear water lakes were found primarily in the interior. Glacial lakes were more evenly distributed but concentrated in watersheds with higher elevations (Fig. 1). These clarity categories were associated with significant differences in limnological conditions and productivity. Euphotic depth, alkalinity, and pH were highest in clear lakes. ANOVAs examining lake

Atlas et al.
Models of sockeye lake rearing capacity biophysical conditions across the three clarity groups revealed that stained lakes had the most acidic water (mean pH = 6.02), compared to less acidic glacial lakes (6.59), and neutral clear lakes (6.95) (p < 0.0001). These physical differences were associated with differences in chlorophyll concentrations among lake clarity groups, with clear lakes having significantly higher concentrations of chlorophyll (p = 0.0046) (Fig. 2). This variation in lake biophysical conditions and regional variations in landscape, climate, and lake clarity were associated with significant differences in mean lake PR. Clear water lakes had the highest mean PR (124.87 mg C m −2 d −1 ), while both glacial (39.22 mg C m −2 d −1 ), and stained lakes (63.87 mg C m −2 d −1 ) had significantly lower mean PR (p < 0.001) (Fig. 3).
Landscape vs. PR relationships AICc supported a range of possible model structures, with 19 different models falling within 4 delta units of the top model (Table 2). Among the landscape and lake-level variables considered, lake-to-watershed ratio (nine models) and perimeter ratio (eight models) appeared in the highest number of models receiving support. Mean watershed slope, % glacial cover, annual precipitation, and the interaction between temperature and water clarity category did not appear in any of the top models suggesting that they explain very little additional variation in the lake PR data. However, AICc ranking do not necessarily reflect the models best suited for prediction. Predictive performance varied among the models selected by AICc, and we further evaluated the predictive performance for the four models with the lowest RMSE using leave-one-out analysis (Supporting Information Fig. S3, Table S2). Table 2. Top ranking linear mixed-effects models of lake PR from AICc model selection, ranked by lowest RMSE. Interactions are denoted by colon. All variables included in interactions are also included as main effects. Water refers to the categorical water clarity group. All models included a random effect of watershed to account for multiple lakes within some watersheds. See Table 1 for definitions of variables. Rankings for models used in prediction are in bold text.

Rank
Model   Water clarity category was included in all the top models ranked by both AICc and RMSE (Table 2). Among the four models with the best predictive performance, the interaction between clarity category and latitude appeared in three of the models, with coefficient estimates (Table 3) suggesting a trend toward lower productivity in higher latitude glacial lakes (−0.52; SE = 0.18), and higher productivity among higher latitude lakes in the stained category (0.17; SE = 0.20). The proportion of the watershed occupied by the focal lake (lake-to-watershed proportion) appeared in three of the best models and had a positive coefficient estimate (0.59, SE = 0.29), indicating higher productivity in lakes which comprise a greater proportion of their watershed area. Perimeter-to-area ratio also appeared in three of the best predictive models, and had a negative effect on lake PR (−0.30, SE = 0.13), such that lakes with longer shorelines relative to their surface area were less productive in the pelagic zone. Upstream lake area appeared in two of the best models and was associated with lower PR values (−0.28; SE = 0.13), and percent tree cover appeared in our fourth best models with higher tree associated with lower PR (−0.29; SE = 0.11).
Watershed level random effects reflect intercept values estimated for each watershed in our sample (n = 49). The degree of variation in these random intercepts reflects part of the difference between the PR in a given watershed relative to the average predicted PR in a watershed with the same clarity category and catchment characteristics. Across the four best predictive models, random intercept values ranged from 5.19 mg C m −2 d −1 in Kitkiata lake to −4.27 in Bloomfield lake (Supporting Information Table S4).

Predicted lake PR
Observed lake PR values ranged from 3.1 to 289 mg C m −2 d −1 reflecting broad variation in lake productivity both within and across lake clarity groups (Supporting Information Table S1). Leave-one-out analysis revealed that across all lakes mean predicted PR fell within 23 mg C m −2 d −1 of the observed value. Given differences in the representation of the three clarity groups in the data set, the uncertainty of model fits and their predictive performance varied by clarity group. In absolute terms, the average difference between predicted and observed PR in clear lakes fell within 33 mg C m −2 d −1 , while stained and glacial lakes had a mean difference in observed and predicted PR of 14 mg C m −2 d −1 and 34 mg C m −2 d −1 , respectively (Supporting Information Fig. S3, Table S2). However, given differences in mean PR among clarity groups, we estimated the relative difference between observed and predicted values by standardizing these differences by the mean PR for a given clarity group. This revealed similar levels of uncertainty in predictions for stained (22%) and clear lakes (26%), and much higher uncertainty for glacial lakes (87%). In particular, the models produced highly uncertain predictions of PR for Owikeno and Kitlope lakes when they were excluded from the sample, suggesting that these lakes exerted considerable influence on the model fit for glacial lakes, likely because they had PR and covariate values that fell outside of the range observed among the sampled lakes.
Across the four models with the strongest predictive performance, predicted lake PR ranged from 3.93 to 165.40 mg C m −2 d −1 , and on average clear lakes were predicted to be almost twice as productive as stained lakes and more than 14 times as productive as glacial lakes. Given the distribution of lake clarity across the NCC-with a high proportion of stained lakes in low elevation coastal watersheds-predicted lake PR showed spatially coherent patterns of productivity across the region (Fig. 4). Among the stained lakes, the highest predicted PR values were in the north in Haida Gwaii and in watersheds around the lower Skeena River, while the lowest predicted PR values were in the more southerly, low elevation watersheds in the Hecate Lowlands. Lakes with the highest predicted PR 2), all clear water lakes. Kimsquit is a highly mountainous drainage in a coastal fjord and has virtually no drainage area above the lake, thus it was a major outlier for the size of the lake relative to the upstream watershed, and predictions for Kimsquit lake that included lake-towatershed ratio produced unrealistically high estimates of PR (Map- Fig. 4, Predictions-Supporting Information Table S5). Among the lakes with the lowest predicted PR, Lower and Upper Kluatantan (7.14; 95% CI = 3.6-14.2) and (8.22; 95% CI = 4.2-16.1) and Oweegee (4.41; 95% CI = 2.3-8.3) were all glacial lakes.
In practical terms, the uncertainty associated with these PR predictions translates into uncertainty for managers seeking to understand the carrying capacity of sockeye rearing lakes. For example, in a stained lake of median size (250 ha) and mean PR for its clarity category, the predicted spawner abundance at carrying capacity (S MAX ) would be 5464 spawners (95% CI = 3798-7972). For glacial and clear lakes of the same size, predicted S MAX in lakes of average productivity would be 3355 (95% CI = 1719-6519) and 10,682 spawners (95% CI = 6758-17,050), respectively. Despite the uncertainty associated with these capacity predictions, they are a major step forward in such as datalimited landscape.  Table S1 and mean predictions and uncertainty for PR in unsampled lakes are presented in Supporting Information Table S5.

Atlas et al.
Models of sockeye lake rearing capacity

Discussion
Canada's Wild Salmon Policy (WSP) was established in 2005 with the goal of protecting wild salmon for the benefit of Canadians in perpetuity. Among the goals of the WSP is the establishment of conservation benchmarks for management and recovery. On the remote NCC of British Columbia, access to many sockeye rearing lakes is difficult, posing significant challenges for population and lake monitoring. As a result, almost 70% of sockeye salmon populations in the region are currently lacking sufficient time series of spawner abundance to evaluate stock status (Pacific Salmon Foundation 2018). To better understand the productivity and juvenile rearing capacity of data-limited sockeye lakes, DFO has conducted rotational limnological sampling in many sockeye rearing lakes across the NCC (Shortreed et al. 1998(Shortreed et al. , 2007 since the 1970s. These monitoring efforts have provided key insights into ecosystem conditions in sockeye rearing lakes, as well as estimates of sockeye carrying capacity that have been used as benchmarks for evaluating current conservation status. However, to date sampling has occurred in fewer than 65 of the 157 sockeye rearing lakes on the NCC. Through our landscape modeling, we generated predictions of PRs for the 96 previously unsampled lakes known to support rearing sockeye salmon. These estimates of PR were then converted to predictions of sockeye lake carrying capacity using the PR model . While management should be grounded in empirical observations of habitat conditions and population size, predictions from our landscape model yielded estimates of carrying capacity with comparable management uncertainty to estimates derived from time series of spawners and recruits (Grant et al. 2011;Connors et al. 2019). Thus, predictions of PR from these 96 unsampled lakes are valuable starting points for the development of conservation benchmarks and can serve as biologically grounded priors for stockrecruit based estimates of management targets.
Our work revealed strong, regionally coherent patterns of sockeye lake productivity across the NCC. Consistent with previous studies, variation in lake biophysical conditions and PR were closely associated with water clarity, reflecting variable light penetration, euphotic volumes, and possibly relative nutrient bioavailability Shortreed et al. 2000Shortreed et al. , 2007, and was further informed in our study by landscape variables at lake and watershed scales. Given the degree to which variation in lake PR was effectively predicted by water clarity and a suite of landscape variables, we generated predictions for 96 sockeye lakes where empirical estimates of lake PR have not previously been made.
Categories of water clarity (clear, glacial, stained) have long been known as an important correlate of lake productivity (Stockner and MacIsaac 1996;Edmundson and Mazumder 2001;Shortreed et al. 2007). Clarity may influence lake productivity through several pathways, including physical limitation via reduced light penetration in humic stained and glacial lakes (Lloyd et al. 1987;Xenopoulos et al. 2003), dampening or amplifying seasonal stratification with effects on lake mixing and nutrient limitation (Stockner and Shortreed 1989;Fee et al. 1996), and mediating the bioavailability of phosphorus and other limiting nutrients (Jackson and Hecky 1980;Edmundson and Carlson 1998;Maranger and Pullin 2002). In the NCC region, the distribution of sockeye nursery lake clarity types is nonrandom, with stained lakes largely concentrated in coastal watersheds where high winter rainfall and peak discharge occur between January and March. Clear lakes typically occupy interior watersheds with snowmeltdominated hydrology and peak discharge during late spring (Stockner andShortreed 1985, 1989). Lacking seasonal ice cover, coastal lakes are mostly monomictic, while interior lakes with continental climates and winter ice cover tend to exhibit dimictic stratification (Stockner 1987;Stockner and Shortreed 1989). Most NCC coastal and interior lakes are oligotrophic, but stained coastal lakes are particularly unproductive. High fall and winter discharge delivers nutrients to coastal lakes during a period when lakes are well mixed and light availability is low, yielding low rates of biological production and limited nutrient retention within lake food webs. Interior lakes by contrast receive peak water and nutrient inputs during the late-spring and early-summer when lakes are undergoing stratification, producing higher rates of nutrient uptake and retention, and higher primary production (Stockner and Shortreed 1985). Glacial lakes are distributed across the longitudinal extent of the NCC. Given their typically mountainous drainages and the contribution of glacial meltwater to their hydrology, glacial lakes typically receive high inputs of water and suspended sediment during spring and summer runoff season, driving physical limitation of biological productivity stemming from high turbidities, shallow euphotic depths, and cold water temperatures (Lloyd et al. 1987;Stockner and MacIsaac 1996;Edmundson and Carlson 1998).
A number of mechanistic drivers may explain the interaction between water clarity and latitude in our study systems, since glacial lake PR was lower at high latitudes, but no similar gradient was apparent for clear or stained lakes. In glacial systems, light is generally less available, and its relative influence on overall seasonal primary production is greater given depthlimited euphotic volumes (Shortreed et al. 2001). Glacial lakes therefore exhibit lower overall seasonal-mean PR than clear water systems (Barouillet et al. 2019). If northern glacial lakes with greater glacial extent and more meltwater inputs have higher turbidity than southern watersheds, elevated turbidity in higher latitude lakes may contribute to the negative relationship between glacial lake PR and latitude. Furthermore, a higher proportion of phosphorus in glacial lakes is nonbiologically available, and the bioavailability of phosphorus may decline in the presence of higher glacial turbidity (Edmundson and Carlson 1998), reinforcing low productivity in lakes with high turbidity. Unfortunately, data are lacking to Atlas et al. Models of sockeye lake rearing capacity evaluate the relationship between turbidity and latitude. Finally, day length shows predictable variation across latitudes, producing shorter growing seasons for northerly lakes. If light penetration is limited in glacial systems, constraining total potential growing season production, glacial lake seasonal mean PR may be inherently more sensitive to changes in light availability than clear or humic stained systems. These linkages between turbidity, light availability, and primary production all provide plausible explanations for the observed pattern of declining glacial lake productivity at higher latitudes. While lake water clarity and latitude explained a large proportion of the variance in lake PR, the geomorphic and climatic attributes of lakes and their watersheds also played an important role in explaining lake productivity. For example, the estimated effect of lake-to-watershed ratio suggests higher productivity in lakes that occupy a larger proportion of their watershed, and the negative effect of perimeter-to-area ratio indicated lower productivity in lakes with more complex shorelines. The size of a lake relative to the upstream catchment has been shown to play an important role in mediating nutrient availability (Prepas et al. 2001), since lakes occupying a large proportion of their watershed area have lower flushing rates and retain water and nutrients longer, permitting autotrophic attenuation (Kratz et al. 1997). Lakes with complex shorelines have more shallow-water littoral habitat supporting benthic algae and macrophytes to comprise a larger proportion of lake-wide primary production and nutrient uptake, potentially limiting the availability of nutrients to limnetic phytoplankton (Jeppesen et al. 1998;Vadeboncouer et al. 2002). In highly oligotrophic systems such as those found in coastal British Columbia, increased primary production and nutrient uptake in the littoral zone may further reduce rates of limnetic production measured by lake PR, particularly given low rates of coupling between littoral and limnetic zones (France 1995). However, lake morphometry, and in particular the presence of more shallow water habitat has been linked to higher fish yields in some humic-stained boreal lakes (Seekel et al. 2018). This source of production is not captured in pelagic-focused modeling (e.g., PR Model- Hume et al. 1996;Shortreed et al. 2000); however, the contribution of littoral primary production and prey such as aquatic insects to juvenile sockeye salmon diets may be considerable in some lakes (Narver 1970;Hume et al. 1996;Richardson et al. 2017). Capacity estimates in lakes with complex shorelines and large areas of shallow water habitat may therefore underestimate total primary production, while still capturing the most important photosynthetic pathway (limnetic) for juvenile sockeye rearing capacity.
Predictions of lake PR serve as meaningful approximations of lake productivity and its influence on secondary production for planktivores; however, lake PR may be more representative of total lake-wide production in some systems than others. Microbial pathways often contribute substantially to production in lake ecosystems (Porter et al. 1988;Weisse and MacIsaac 2000), particularly in highly stained lakes where high rates of microbial production can produce net heterotrophy (Jansson et al. 2000;Ask et al. 2009). Food web structure and community composition often differentiate strongly between stained and clear lakes in the coastal British Columbia and Alaska (Stockner 1987;Stockner and Shortreed 1989;Koenings et al. 1990). With strong nutrient limitation and low rates of autotrophic production, stained lakes typically have higher concentrations of picoplankton, and a greater dependence on microbial pathways for basal production (Stockner 1987;Stockner and Shortreed 1989). Given their small size, the high abundance of picoplankton and bacteria in stained lake food webs adds additional trophic levels between primary production and planktivorous fish, making energy transfers to higher trophic levels much less efficient (Stockner and Shortreed 1989). Unproductive stained and glacial lakes often lack Daphnia and other large-bodied cladocerans, with limnetic grazer communities dominated by rotifers and other smaller-bodied zooplankton (Stockner and Shortreed 1989;Koenings et al. 1990). Small zooplankton may serve as an energy sink if they are too small to be consumed by planktivorous fish (O'Neill and Hyatt 1987;Stockner and Shortreed 1989). These differences in food web structure ultimately reduce the amount of energy available to higher trophic levels in the limnetic food webs of stained lakes. Community composition, food web structure, and the magnitude of heterotrophic energy flows therefore represent an important and currently unquantified contributor to lake energy budgets, constituting a majority of lake-wide production in some instances (Nürnberg and Shaw 1998). Changes in nutrient availability (Weisse and MacIsaac 2000) and DOM inputs (Jansson et al. 2000) can act to modulate the importance of microbial pathways in lake food webs. Given the multiple interacting trophic pathways supporting coastal lake food webs, models accounting for heterotrophic production would likely improve the accuracy of estimated sockeye lake rearing capacity in heavily stained coastal lakes.
Despite the advances offered by landscape-scale predictive models, ongoing limnological assessments are a vital component of managing lake-dependent sockeye populations. Lake photosynthesis and total food web productivity can exhibit both directional change and interseasonal to interannual variability (Fee 1980;D. T. Selbie unpubl.) in response to climate and hydrologically mediated variability in the delivery of nutrients and organic material (Jansson et al. 2000), and changes in the biogenic delivery of nutrients via spawning salmon (Stockner and MacIsaac 1996;Schindler et al. 2005;Chen et al. 2011). Climate change may therefore drive changes in nutrient dynamics and lake productivity in unexpected ways (Adrian et al. 2009) as the hydrology of many systems transitions from snowmelt to rain-dominated (Klos et al. 2014), and ongoing declines in the survival of salmon in the ocean reduce the delivery of salmon derived nutrients to coastal watersheds (Larkin and Slaney 1997). Furthermore, glacier retreat may lead to changes in lake productivity in the immediate future (Bliss et al. 2014). In light of the dynamic nature of sockeye rearing lakes, continued monitoring of lake food web productivity and structure will provide necessary insight into the physical and biological conditions that drive the freshwater population dynamics of sockeye salmon.
Landscape ecology perspectives and approaches are increasingly being applied to aquatic systems, with researchers and conservation practitioners seeking to understand the influence of spatial patterns, landscape context, and linkages between adjacent habitat patches for aquatic ecosystem structure and function (Wiens 2002). Landscape-scale models have been proposed for classification and management of aquatic ecosystems in regions where monitoring and managing individual lakes or stream ecosystems may be infeasible (Soranno et al. 2010). These approaches are tailored to the needs of resource managers attempting to inform management across broad spatial scales in data-limited landscapes. Given the remote nature of NCC lakes and the cost and logistical challenges associated with monitoring, generating predictive models of ecosystem productivity or function can inform conservation and management, and guide investments in more intensive monitoring. Model outputs will therefore provide interim inputs for PR model estimates of lake juvenile rearing capacity, which can be improved through future sampling of lake trophic structure and productivity.
Our results demonstrate the close links between lake and landscape attributes and lake primary productivity, yielding predictions of lake PR for remote watersheds across the NCC which have previously been unsampled. Given previous research linking lake PR to sockeye salmon juvenile rearing capacity , our findings and the resulting predictions of lake PR will provide a starting point for evaluating the productive potential of sockeye rearing lakes across the NCC. However, these predictions should not be viewed as a substitute for robust limnological and population monitoring, which are essential for precautionary management of fisheries and detection of environmental changes. Efforts are currently underway to assess the status of data poor sockeye populations in British Columbia and understanding the links between lake and landscape characteristics and lake PR will serve as a critical stepping stone toward evaluating conservation status and developing data-driven management approaches for sockeye populations with limited data.