Here, we generate lake‐specific predictions of walleye recruitment and largemouth bass relative abundance in Wisconsin lakes under a future climate scenario. We use a mechanistic, thermodynamic simulation model to predict daily, depth‐specific water temperatures for over 2100 Wisconsin lakes in contemporary (1989–2014), mid‐century (2040–2064), and late‐century (2065–2089) climatic conditions. Mechanistic, thermodynamic models of water temperature can capture heterogeneity in lake responses to climate conditions, but to date, such models have been applied to either single lakes (e.g., Hadley et al ., 2013 ), or to ‘model’ lakes across the landscape (e.g., Fang & Stefan, 2009 ). We link contemporary modeled water temperatures and other lake characteristics to the reproductive success of walleye and relative abundance of largemouth bass, and project changes in both species in two future periods. In doing so, we provide a framework for adaptation to climate change on a scale and resolution relevant to resource management.

Air temperature is often used as a proxy for water temperature when projecting the effects of climate change on fish species (e.g., Chu et al ., 2005 ; Sharma et al ., 2007 ; Hein et al ., 2011 ; Van Zuiden et al ., 2016 ). Water temperature and air temperature are strongly correlated across large latitudinal gradients such as those encompassed by, for example, Canada (Chu et al ., 2005 ; Sharma et al ., 2007 ; Van Zuiden et al ., 2016 ) or Sweden (Hein et al ., 2011 ). Air temperature can be a good proxy for an integrated metric of lake temperature such as mean summer surface temperature (Sharma et al ., 2008 ). Using air temperature as a proxy assumes that all lakes in a region respond similarly to similar climate warming. However, lakes respond differently to climate depending on their size, clarity, depth, and sheltering from wind (Fang & Stefan, 2009 ; Read et al ., 2014 ; O'Reilly et al ., 2015 ). Consequently, air temperature alone may not capture differences among lakes in seasonal and vertical temperature patterns, factors known to be important in structuring lake fish communities (Magnuson et al ., 1979 ). Thus, using air temperature to infer current and future lake temperatures may not adequately characterize relevant components of thermal habitat for fish and how they might change with climate warming. Accounting for local‐scale heterogeneity in environmental conditions is critical for models of species distribution and population processes such as growth and reproduction (Dormann et al ., 2012 ). Approaches that account for heterogeneous responses of local habitat to climate can improve the accuracy of predicted species responses to climate change by allowing for local‐scale variability that is not possible when predictions are made at coarse spatial resolutions. Such approaches can identify refugia where species are expected to persist in spite of unfavorable regional conditions, and these refugia are critical for adaptation to climate change (Suggitt et al ., 2011 ; Keppel & Wardell‐Johnson, 2012 ).

Fish communities in Wisconsin lakes have changed in the past several decades; the reproductive success of walleye ( Sander vitreus ; a coolwater species) has declined, while the abundance of largemouth bass ( Micropterus salmoides ; a warmwater species) has increased (Hansen et al ., 2015b , c ). Both species are frequently targeted by anglers, although walleye are generally the preferred species (Mcclanahan & Hansen, 2005 ; Holsman et al ., 2016 ). Walleye are more likely to successfully reproduce in lakes with cooler water, but this effect of water temperature varies with lake size and other characteristics (Hansen et al ., 2015a ), suggesting that some lakes may be more resilient than others to future changes in climate. Identifying lakes that are likely to support natural walleye recruitment now and in the future could allow resource managers to adapt to climate change by guiding the allocation of resources to lakes most likely to continue to support desired species, and through outreach with stakeholders to set reasonable expectations for fish communities.

Fish communities in north temperate lakes are composed of species with distinct temperature preferences (i.e., warm‐, cool‐, and coldwater species), and all three thermal guilds can coexist in thermally stratified lakes (Magnuson et al ., 1979 ). Many studies that project the effects of climate change on lake fishes have predicted declines in coldwater species (Fang et al ., 2004b ; Mackenzie‐Grieve & Post, 2006 ; Jacobson et al ., 2010 ; Herb et al ., 2014 ) and increases in warmwater species (Lehtonen, 1996 ; Chu et al ., 2005 ; Sharma et al ., 2007 ; Van Zuiden et al ., 2016 ) as climate warms (reviewed in Comte et al ., 2013 ). Coolwater fish thermal habitat frequently is also predicted to increase under climate change (Buisson et al ., 2008 ), especially in north temperate areas such as the Great Lakes region of North America (Fang et al ., 2004a ; Jansen & Hesslein, 2004 ; Chu et al ., 2005 ). However, predictions are not consistent across regions and species (Comte et al ., 2013 ; Van Zuiden et al ., 2016 ). In addition, temperature influences fish through numerous pathways, including growth, metabolism, feeding rates, reproductive timing, and more (Kitchell et al ., 1977 ; Magnuson et al ., 1979 ). Recent work has linked failed reproductive success of a Great Lakes coolwater species to reduced winter severity and shifts in the timing of temperature‐driven spawning cues (Farmer et al ., 2015 ). Therefore, the effects of climate change on coolwater fish species may be complex (Shuter et al ., 2012 ), and may not be captured by changes in preferred thermal habitat alone.

Climate is a primary determinant of species distribution and abundance, and projecting future species distributions under climate change is an important goal for conservation and management (Hannah et al ., 2002 ; Heller & Zavaleta, 2009 ). Freshwater ecosystems are particularly vulnerable to the effects of future climate change (Millenium Ecosystem Assessment, 2005 ). In particular, lakes serve as ‘sentinels’ of climate change by integrating effects through a variety of pathways (Adrian et al ., 2009 ; Schindler, 2009 ; Williamson et al ., 2009 ). Many studies have used present‐day distributions of species correlated with environmental conditions to project climate‐driven extinction risk and changes in distributions of freshwater taxa at a variety of spatial extents and resolutions (reviewed in Comte et al ., 2013 ; Heino et al ., 2009 ). However, predictions of relative abundance or natural reproduction of important species will be more useful for fisheries managers than predictions of presence/absence. Additionally, decision makers operating on local (e.g., state‐level) scales will benefit from predictions made at small (e.g., lake‐specific) spatial resolutions to enable effective decision making about resource allocation and management actions.

Using contemporary fish models for each species, we predicted the probability of walleye recruitment success and high largemouth bass abundance for all lakes in the contemporary period (including lakes never surveyed for walleye or largemouth bass). We also predicted future lake‐specific probabilities of walleye recruitment success and high largemouth bass abundance using projected values of thermal metrics in mid‐ and late century from the three GCMs, while holding nonthermal variables constant. We applied thresholds identified from the contemporary model to classify lakes as having successful/unsuccessful walleye recruitment and high/low abundance of largemouth bass in the two future periods. Individual lakes were classified as those predicted to gain, lose, or maintain each species. Finally, we classified lakes based on the median predicted probabilities of the two species in each time period across the three GCMs as follows: walleye dominant (successful walleye recruitment, low largemouth bass abundance), bass dominant (failed walleye recruitment, high largemouth bass abundance), coexistence (successful walleye recruitment, high largemouth bass abundance), or neither (failed walleye recruitment, low largemouth bass abundance). All results are based on lakes for which all predictor variables for the two fish models were available for each time period ( n = 2148).

Following identification of the final model for each species, we used multiple metrics of model accuracy in classifying fish abundance in the contemporary period. Area under the curve (AUC) ranges from 0.5 for models that are no more accurate than expected by chance to 1.0 for models that are 100% accurate in classifying binary outcomes (Manel et al ., 2001 ). AUC is useful as a measure of accuracy because predicted probabilities need not be converted to binary responses based on an artificial threshold value, and AUC is not biased by species prevalence – or in our case, recruitment success or high largemouth bass abundance – as are other accuracy metrics (Fielding, 1999 ; Manel et al ., 2001 ). We identified probability thresholds separating failed from successful recruitment and low from high largemouth bass abundance; thresholds were selected to maximize sensitivity (true predictions of high abundance or recruitment success) plus specificity (true predictions of low abundance or recruitment failure) using the PresenceAbsence package in r (Freeman & Moisen, 2008 ). As additional measures of accuracy, we also report the proportion of cases in which the model classified lakes correctly in the contemporary period. These accuracies are reported both in terms of overall model accuracy and class‐specific accuracies measured using OOB data in the randomForest package.

We used a variable selection procedure to balance predictive accuracy with model simplicity, adapted from Díaz‐Uriarte & De Andres ( 2006 ) and described in detail in Hansen et al . ( 2015a ). For each species, we fit 25 random forest models that included all potential predictor variables (Tables S1 and S2), and calculated variable importance as the mean of the mean decrease in accuracy associated with each variable (Cutler et al ., 2007 ) across all 25 models. Overall model error (proportion of inaccurate classifications) was also averaged across the 25 models. Variables were then dropped sequentially in reverse order of importance, and 25 random forest models were fit using the remaining variables with each subsequent drop. Overall model error rate was calculated with each subsequent drop. The simplest model (with the fewest predictor variables) with an overall error rate within 1 SE of the minimum error rate across all models was selected as the final model.

We used random forest models to predict walleye recruitment success and largemouth bass relative abundance in Wisconsin lakes, implemented in the randomForest (Liaw & Wiener, 2002 ) package in r version 3.2.1 (R Development Core Team, 2015 ). Random forest models identify interactions and nonlinear relationships that need not be specified a priori (Breiman, 2001 ), and are excellent tools for prediction (Breiman, 2001 ; Prasad et al ., 2006 ; Cutler et al ., 2007 ). Furthermore, random forest models have proven useful for predicting walleye recruitment success in Wisconsin lakes (Hansen et al ., 2015a ). Multiple (1000) classification trees comprise random forest models, and input data are bootstrapped for the construction of each tree. Approximately 2/3 of the full dataset are used for the construction of each tree, with the remaining data used as the validation dataset to assess the out‐of‐bag (OOB) error rates. Model accuracy is measured as average error rate in predicted class membership across all trees in the forest based on OOB data. Only a subset of predictor variables are evaluated for each potential split point in the tree, which avoids model overfitting even with a large number of potential predictor variables (Cutler et al ., 2007 ). The influence of each predictor variable (mean decrease in accuracy) is measured by randomly permutating predictor variable values and calculating the mean reduction in model accuracy across all trees that results from using randomly permutated values in place of true values. We used the stratification feature of the randomForest package to avoid pseudoreplication (Pecl et al ., 2011 ). Because we had multiple years of data for our response variables (recruitment success or largemouth bass relative abundance) in some lakes but our predictors were constant for each lake, we stratified bootstrapped samples by lake to ensure that each tree was constructed using a maximum of one observation from each lake. The predicted response from classification random forest models is the probability of belonging to a specific class (e.g., probability of successful walleye recruitment). We set the terminal node size to five to ensure that multiple lakes were included in each tree.

Lake‐specific predictors of walleye recruitment success and largemouth bass relative abundance included measures of lake morphometry, productivity, and species‐specific thermal habitat metrics (Tables S1 and S2). When multiple predictor values were available for a single lake, mean values for the contemporary period (1989–2014) were used. We did not include any land‐use or land‐cover variables, because they were not identified as important predictors of walleye recruitment success in a previous analysis (Hansen et al ., 2015a ), and because land‐use projections were beyond the scope of this project. Where possible, we used direct variables in place of indirect (Guisan & Thuiller, 2005 ; Austin & Van Niel, 2011 ). For example, temperature directly affects fish metabolism, growth, and survival, while latitude serves as a less specific proxy. If two predictor variables were highly correlated (| r | ≥ 0.8), we excluded one of the two from the initial base model, when possible retaining the more ecologically relevant variable (Tables S1 and S2).

We converted quantitative metrics of relative abundance of each species to binary data representing successful or failed recruitment for walleye and high or low abundance for largemouth bass. We avoided quantitative measures of relative abundance for several reasons. Fisheries catch rates are highly variable from year to year because of both variability in fish populations and sampling errors (Ricker, 1975 ; Hilborn & Walters, 1992 ); as a result, continuous abundance values are difficult to predict with any degree of certainty. Still, electrofishing provides a useful index of relative abundance (Schoenebeck & Hansen, 2005 ), particularly on the coarse scale of ‘high’ vs. ‘low’ abundance that is of interest here. We translated continuous catch rates into binary classes using expert opinion of fisheries managers (Fernandes et al ., 2010 ). Walleye recruitment was defined as successful if age‐0 walleye catch rates in fall electrofishing surveys were greater than 6.2 fish per km (10 per mile; Hansen et al ., 2015b ), and low otherwise. Largemouth bass relative abundance was classified within seasons (spring and fall) as ‘high’ if greater than the season‐specific global median and ‘low’ otherwise, and combined into one dataset. Model results for both species were insensitive to changes ±30% in these thresholds (data not shown).

Fish data were collected by the WDNR and the Great Lakes Indian Fish and Wildlife Commission as a part of their standardized monitoring programs. Walleye reproductive success was measured as successful recruitment to age 0 from naturally reproduced fish. This walleye recruitment metric was indexed by fall electrofishing catch of age‐0 walleye per surveyed kilometer collected in nighttime fall surveys using a 230 V AC electrofishing boat. Age‐0 walleye were identified either by length–frequency analysis or using scale annuli following WDNR protocols (Serns, 1982 ; Hansen et al ., 2004 ). We used lake survey data from 1989–2014, excluding surveys from lake‐years where age‐0 walleye were stocked prior to the survey. To ensure that effort was sufficient to index whole‐lake abundances of age‐0 walleye, we restricted our analysis to surveys where at least 70% or 16.1 km of shoreline was sampled, as per WDNR age‐0 walleye sampling protocols. Largemouth bass relative abundance was indexed from 1995–2014 using catch per kilometer of largemouth bass over 20.3 cm (8 inches) collected in both fall and spring using a 230 V AC electrofishing boat. In some largemouth bass surveys, time instead of distance shocked was recorded as the measure of survey effort. In these cases, hours were converted to miles based on the empirical relationship observed in the data (Effort miles = 2.035 × Effort hours + 0.044, R 2 = 0.84, P < 0.0001, n = 9043; J. F. Hansen, personal communication ). For both walleye and largemouth bass, we excluded surveys in which water temperatures were <10 °C or >21 °C (following WDNR survey protocols) and those in which survey reliability was low because of adverse conditions identified by the survey operator. Initially, we analyzed 2608 walleye recruitment values from 470 lakes and 3350 largemouth bass relative abundance values from 1126 lakes, although final sample size depended on the results of the model selection process (described below) because some predictor variables were not available for all lakes.

In this study, it was critical to predict more than simple presence and absence in order to be relevant to management decision making. Largemouth bass are present at some level in most Wisconsin lakes, but whether population densities are high or low is most relevant to fisheries management and fish community dynamics. Largemouth bass population density is commonly indexed as relative abundance measured as electrofishing catch per effort (Schoenebeck & Hansen, 2005 ; Reynolds & Kolz, 2012 ). For walleye, adult abundance can be a misleading response variable, as many lakes contain walleye populations supported entirely by stocking. In general, populations supported by natural recruitment reach higher densities than those supported by stocking (Nate et al ., 2000 ), and management strategies for maintaining naturally reproducing populations are likely to be quite different than those relevant for stocked lakes. Additionally, natural reproduction is declining in Wisconsin lakes (Hansen et al ., 2015b ). Thus, we used reproductive success as the relevant metric for walleye relative abundance.

Daily lake temperature profiles from the four meteorological driver datasets (NLDAS in contemporary period and CM2.0, ECHAM5, and GENMOM for future periods) were used to derive species‐specific annual temperature metrics hypothesized to be biologically relevant based on published temperature preferences or correlations with each fish species (Tables S1 and S2). We calculated a number of metrics describing the average or cumulative surface water temperature of a lake, including water temperature degree days with a base temperature of 5°C (Tables S1 and S2; Chezik et al ., 2013 ). Multiple temperature ranges and cutoffs were used to account for various definitions of optimal thermal habitat that apply to various life stages of each species (e.g., Coutant, 1977 ; Wismer & Christie, 1987 ). We also attempted to quantify relevant thermally driven processes for fish, such as temperature‐driven spawning cues or winter duration for egg production (Shuter et al ., 2012 ; Farmer et al ., 2015 ). All thermal metrics were calculated in r version 3.2 (R Development Core Team, 2015 ) using the rLakeAnalyzer (Read et al ., 2011 ; Winslow et al ., 2015 ), insol (Corripio, 2014 ), and plyr (Wickham, 2011 ) r packages. Code for calculating each thermal metric is available online ( https://github.com/USGS-R/mda.lakes/ ). Lake‐specific mean values of each thermal metric were calculated for each of the three study periods – contemporary (1989–2014), mid‐century (2040–2064), and late‐century (2065–2089) – and used as inputs for statistical models of fish abundance. All modeled temperature metrics for each lake and time period are available online (Hansen et al ., 2016 ).

We adjusted the downscaled GCM meteorological driver data to eliminate bias, because the original data source (Hostetler et al ., 2011 ) was not de‐biased to locally observed meteorology. To do this, we compared downscaled meteorology with NLDAS data for the years where datasets overlapped in the contemporary period (1980–1999). For our study region, air temperature, shortwave radiation, and wind speed were adjusted using different methods depending on the nature of the bias during the overlapping contemporary period. Air temperature was adjusted by creating a linear model between daily GCM and NLDAS average air temperature and applying the modeled parameters to future air temperatures. Shortwave radiation was adjusted by subtracting the offset between mean shortwave for the GCM and NLDAS for each lake. Wind speed was corrected using a multiplicative correction based on the quotient of mean wind speed for each GCM and NLDAS during the contemporary overlap. Each GCM projects future conditions over a different time period; CM2.0 included 2040–2069, ECHAM5 included 2020–2099, and GENMOM included 2020–2089. As with NLDAS, these data were processed using the geoknife r package for each lake. Temporal averaging of meteorological drivers was not used because the downscaled climate models provided only daily output. Lakes for which temperature simulations failed for any GCM were excluded from further analysis.

Future lake temperatures were projected using the same methodology as the hind‐casted simulations, except that three data sources were used for meteorological driver data in place of NLDAS data. Dynamically downscaled climate drivers produced by Hostetler et al . ( 2011 ) were used from the CM2.0, ECHAM5, and GENMOM (Adler et al ., 2011 ) global climate models (GCMs). These three GGMs encompass a wide range of air temperature projections; air temperature projections under GENMOM are low relative to other GCMs, those under ECHAM5 are near the median (Herb et al ., 2014 ), and those under CM2.0 are relatively high (Hostetler et al ., 2011 ). Details for the downscaling procedure used can be found in Hostetler et al . ( 2011 ). Each GCM was driven by the A2 future scenario (Nakicenovic et al ., 2000 , IPCC, 2007 ). The A2 scenario is a pessimistic estimate of human response to greenhouse gas emissions, and predicts the second highest level of warming by the end of the 21st century across all six IPCC scenarios. A2 is similar to the updated upper bound scenario (RCP8.5) developed for the IPCC 5th Assessment Report (Riahi et al ., 2007 ; Hostetler et al ., 2011 ).

For validation of modeled lake temperatures in the contemporary period, we compared model outputs to a large‐scale database of water temperature observations provided by the Wisconsin Department of Natural Resources (WDNR) and the North Temperate Long‐term Ecological Research site ( http://lter.limnology.wisc.edu ) collected as part of their respective standard monitoring programs. The water temperature database consists of 274 348 water temperature observations across 1152 lakes sampled from 1967 to 2013, although data not overlapping our simulation period were ignored. For validation, observations were paired with modeled estimates for the same lake, date, and depth, and root‐mean‐square error (RMSE) was calculated. For further understanding of the depth dependence of model uncertainty, volumetrically averaged epilimnetic and hypolimnetic temperatures were calculated from observed and modeled profiles using rLakeAnalyzer (Read et al ., 2011 ; Winslow et al ., 2015 ), and those layer‐average errors were also summarized using overall RMSE.

Similar to Read et al . ( 2014 ), lake temperatures were ‘hind‐casted’ for the contemporary time period (1989–2014) using gridded meteorological drivers provided by the North American Land Data Assimilation System (NLDAS, Mitchell et al ., 2004 ). Time series inputs of downwelling long‐ and shortwave radiation, air temperature, wind speed, relative humidity, and precipitation drive the lake model. The boundary of each lake was used to create an area‐weighted average of each driver from the NLDAS gridded data; data were processed using the geoknife r package (Read et al ., 2015 ). Hourly NLDAS data were converted to daily means for each lake, except wind speed data that were power‐averaged to better represent the average energy in wind, . Wind speed values during the contemporary period were adjusted to remove the influence of a step change in the data that began in 2001, which likely resulted from a change in NLDAS from retrospective to real‐time data assimilation (Mitchell et al ., 2004 ). We adjusted wind speeds by multiplying raw wind speed data by 0.921 for all wind speed data after December 31, 2001, based on the ratio between the scale parameters for Weibull distributions fit to pre‐ and post‐2001 NLDAS wind data. The correction was applied post‐2001 as comparisons with local meteorological station data showed pre‐2001 data closely matched the observed wind speed distribution.

We used the lake modeling methodology documented in Read et al . ( 2014 ), with two notable changes. First, Read et al . ( 2014 ) used an empirical model to define ice‐on and ice‐off dates, and water temperatures were simulated only for the ice‐free season for each lake. Here, we used an improved snow/ice module in GLM (implemented in GLM version 2.1.0) to model ice cover, define ice‐free periods, and run continuous simulations during the study periods. Second, water clarity data used to parameterize the light extinction properties of each lake were selected hierarchically in this study to use the best available data for each lake. When available, we used the mean of all lake‐specific in situ Secchi depths (1552 lakes). When in situ measurements of Secchi depth were not available, we used the mean of all lake‐specific satellite‐estimated Secchi depths (methods described in Torbick et al ., 2013 ; 579 lakes). In contrast, Read et al . ( 2014 ) used a lake‐specific average of all observed Secchi depth values regardless of data source ( in situ or satellite). All other lake temperature modeling routines used in this study followed Read et al . ( 2014 ).

We hind‐casted and forecasted water temperatures of all lakes in Wisconsin, US, with a history of any fisheries management activity (e.g., a fisheries survey or stocking event) for which lake‐specific data were available ( n = 2148; data needs described below). Daily lake temperature profiles for all study lakes were estimated using the General Lake Model (GLM) version 2.1.8 (Hipsey et al ., 2013 ). This model is an open‐source, mechanistic lake model based on algorithms that were developed for earlier one‐dimensional mechanistic lake models (Hamilton & Schladow, 1997 ). The model is parameterized using lake‐specific properties including water clarity, morphometry, and wind sheltering (see Hipsey et al ., 2013 for a list of parameters and their meanings). To reduce model instability, the model was run at an hourly time step with daily mean meteorological drivers. Only daily outputs were used in results. Lakes for which simulations failed in any year were excluded from further analysis.

Number of lakes classified by predicted species dominance under contemporary (1989–2014), mid‐century future (2040–2064), and late‐century future (2065–2089) conditions. Lakes were classified based on the median predicted walleye recruitment success and largemouth bass relative abundance across three global circulation models for mid‐century climate and two global circulation models for late‐century climate. A total of 2148 lakes are represented here; these are the lakes for which all variables used in the predictive models of both species were available for each time period. Classes are defined as: walleye dominant (blue; walleye recruitment success and low largemouth bass relative abundance), coexistence (purple; walleye recruitment success and high largemouth bass relative abundance), largemouth bass dominant (orange; walleye recruitment failure and high largemouth bass relative abundance), or neither (gray; walleye recruitment failure and low largemouth bass relative abundance). Colored lines show movement among lake classes, and line width is proportional to the number of lakes moving between each class.

Future climate change was predicted to shift fish species dominance in many Wisconsin lakes based on median projected probabilities of supporting each species across GCMs (Fig. 3 ). Of the 2148 lakes for which predictions of both species were possible for each time period, 1236 (58%) were predicted to be dominated by largemouth bass (i.e., relative abundance of largemouth bass predicted to be high, and walleye recruitment predicted to fail) under contemporary conditions. The percentage of lakes dominated by largemouth bass was predicted to increase to 1857 (86%) and 1961 lakes (91%) under mid‐ and late‐century climate conditions, respectively. Conversely, the number of lakes dominated by walleye was predicted to decline from 184 lakes (8.6%) in contemporary conditions to only 25 (1.2%) and 17 lakes (0.8%) in mid‐ and late‐century conditions, respectively (based on median predicted probabilities across GCMs; Fig. 3 ). The majority (58%) of present‐day walleye lakes were predicted to transition into bass‐dominated lakes by mid‐century, although 21% of present‐day walleye lakes were predicted to support both species in the future. Interestingly, the number of lakes predicted to support both species increased slightly over time, although the identity of these lakes changed (Fig. 3 ). Most lakes where coexistence was predicted for contemporary conditions were predicted to become dominated by bass in the future, but the lakes transitioning from walleye‐dominant to coexistence lakes more than offset this loss. Lake‐specific median predicted probabilities of supporting each species under the three time periods and their associated fish dominance categorizations are available in Table S3, and as a shapefile online (Hansen et al ., 2016 ).

The probability of successful walleye recruitment was predicted to decrease and the probability of high largemouth bass abundance was predicted to increase in most lakes based on predicted mid‐ and late‐21st century climate conditions (Fig. S3). As a result, from 74 to 168 Wisconsin lakes that currently are capable of supporting successful walleye recruitment were predicted to lose this capacity (depending on time period and GCM; Table 3 ). Furthermore, from 335 to 771 lakes were predicted to gain the capacity to support high relative abundances of largemouth bass (again depending on time period and GCM). For both species, predicted changes were smallest under the GENMOM model in mid‐century conditions, while all other models and time periods predicted similar outcomes. Importantly, even under late‐century conditions, from 72 to 85 lakes were predicted to support successful walleye recruitment (Table 3 ).

Degree days (the only thermal variable in the final fish models) were projected to increase in Wisconsin lakes by the mid‐ and late‐21st century, although the magnitude of change varied among lakes, climate models, and time periods (Fig. 2 ). Mid‐century predictions ranged from 4% decreases to 17% increases in degree days; late‐century predictions ranged from 2% increases to 31% increases in degree days (Table 2 ). The GENMOM model predicted the lowest degree day increases under future climate conditions; under this model, the range of expected degree days in mid‐century completely overlapped the contemporary range (Fig. 2 ). Higher degree days were predicted in late century compared to mid‐century, although again, GENMOM model predictions were less extreme than the ECHAM5 model (projections from the CM2.0 model were not available for this time period).

Water temperature degree days were correlated with both species, but affected walleye recruitment and largemouth bass relative abundance in opposite ways (Fig. 1 ). Walleye recruitment was most likely in lakes with low degree days, while high largemouth bass abundance was most likely in lakes with high degree days. Responses of both species changed abruptly at a threshold of 2400–2500 degree days. The mirror image relationship between degree days and the two species suggests that lakes that are most suitable for walleye recruitment are least likely to have high largemouth bass abundance; conversely, lakes with high largemouth bass abundance are least likely to support natural walleye recruitment.

Predictions of largemouth bass relative abundance were accurate when validated using independent data (Table 1 ). The AUC value (0.86) indicates high accuracy independent of the threshold used to separate high from low largemouth bass abundance. Using the optimal threshold that maximizes prediction accuracy (0.49), largemouth bass relative abundance was predicted correctly in 79% of lake‐years. Class‐specific accuracies were nearly identical when predicting largemouth bass relative abundance, although high largemouth bass abundance was predicted slightly more accurately than low abundance (Table 1 ). In other words, our model is slightly more likely to generate false negatives (predicting low bass abundance where it is actually high) than false positives (predicting high bass abundance where it is actually low).

Contemporary predictions of walleye recruitment were accurate based on a number of criteria (Table 1 ). The AUC of 0.84 indicates high model accuracy independent of any threshold used to separate successful from failed recruitment (Manel et al ., 2001 ). Other metrics of classification success are contingent upon the threshold used to separate successful from failed recruitment. Using the optimal threshold that maximized prediction accuracy (0.49), walleye recruitment success was predicted accurately in 80% of cases. Recruitment success was predicted more accurately than recruitment failure (Table 1 ). In other words, our model was more likely to generate false negatives (predicting failed recruitment in instances where recruitment was successful) than false positives (predicting recruitment success in instances where recruitment failed).

Modeled water temperatures were accurate when compared to observed data. Comparing all observed water temperatures (all seasons, all depths), the overall root‐mean‐square error (RMSE) was 2.54 °C ( n = 274 348). Volumetrically averaged epilimnion temperatures were more accurate (RMSE = 1.87 °C; n = 23 700). Surface temperatures, defined as temperatures at depths 0–2 m, had slightly lower accuracy (RMSE = 1.94 °C, n = 130 560) than average epilimnetic temperature. Despite using default model parameters across all lakes, our estimated RMSE is only about 1 °C higher than a recent multilake study that used mechanistic thermodynamic simulations calibrated to individual lakes (Fang et al ., 2012 ). Models errors did not have a systematic bias, with an average difference between modeled and observed temperatures of −0.03 °C. Mean degree days in the contemporary period varied between lakes, ranging from 2098 to 3777 °C days, with a median value of 2540 °C days.

Discussion

Climate change undoubtedly will influence the species composition of temperate lakes, with many lakes becoming more suitable for warmwater species and less suitable for cool‐ and coldwater species. Indeed, these changes may already be underway (Farmer et al., 2015; Lynch et al., 2016). However, fish community responses to climate change are complex due to heterogeneous lake responses to climate change and the myriad pathways through which temperature influences fish. We predicted that thermal regimes of Wisconsin lakes will not respond homogeneously to increased temperatures under a high‐end warming scenario. Furthermore, although the relative abundances of two important fish species were influenced strongly by water temperature, the effect of temperature varied among lake types. As a result, we identified resilient lakes that we predicted will continue to support natural walleye reproduction despite a warming climate. These results could help guide the allocation of management resources and tailor specific management actions in locations where they are most likely to succeed.

Our results demonstrate heterogeneity in both lake responses to climate change and the effects of lake temperature on habitat suitability for fish species. Contemporary thermal regimes of Wisconsin lakes differ with lake size, depth, clarity, and wind sheltering (Read et al., 2014), and projected future thermal regimes were predicted to be similarly diverse even within a single GCM and climate scenario (Fig. 2). Furthermore, walleye and largemouth bass populations in diverse lake types did not respond equally to temperature. Although both species were influenced by water temperature degree days, with nearly mirror image functional responses on average, the effect of degree days varied among lakes (Fig. 1). Future warming was predicted to negatively influence walleye recruitment. Only 3.3–4.0% of Wisconsin lakes were projected to support natural walleye recruitment by the late 21st century (Table 3), with most present‐day walleye lakes becoming more suitable for largemouth bass. Conversely, future climate conditions were predicted to positively influence largemouth bass abundance; from 89–95% of Wisconsin lakes were predicted to support high largemouth bass abundance by late 21st century (Table 3).

Some lakes were predicted to be resilient to climate warming; we identified up to 85 lakes predicted to maintain natural walleye recruitment and up to several hundred lakes predicted to maintain low largemouth bass abundance even under an extreme emission scenario (lake names and locations of these resilient lakes are listed in Table S3). These lakes represent local refugia where thermally sensitive species such as walleye can thrive even in regionally unfavorable climate (Ashcroft et al., 2009). These refuge lakes are critical for the persistence and management of species negatively influenced by warming temperatures, but such refugia are only identifiable using approaches that account for small‐scale heterogeneity in the responses of local environments (i.e., lakes) to climate change (Suggitt et al., 2011; Keppel & Wardell‐Johnson, 2012; Keppel et al., 2015). Predictions of future fish community composition based on air temperature alone are unlikely to identify such refugia. The same is true for streams, where approaches that consider the effects of local groundwater inputs identify heterogeneous responses to climate that are not apparent based on coarse‐scale regional air temperature projections alone, with important implications for the projected effects of climate change on fishes (Lyons et al., 2010; Snyder et al., 2015).

Our fish distribution models based on water temperature degree days, lake morphometry, and productivity (conductivity for walleye and Secchi depth for largemouth bass) were accurate based on a number of metrics (Table 1); our model accuracy and AUC values equal or exceed many in the published literature (e.g., Chu et al., 2005; Sharma et al., 2007; Buisson et al., 2008; Hein et al., 2011). Although the lake temperature model predicts entire water column temperatures, in this case, a measure of surface water temperature was the best thermal predictor of both species. We assessed the predictive power of a variety of temperature metrics hypothesized to influence each species, and in doing so, identified surface water degree days (base temperature 5 °C) as the best temperature‐derived predictor of both species. Interestingly, our temperature models predicted increased available thermal habitat for walleye in most lakes when habitat is measured as the proportion of the water column or number of days per year where water temperatures fall between 18–22 °C (results not shown), similar to the increase in coolwater fish habitat predicted by other studies (e.g., Fang et al., 2004a; Chu et al., 2005). By assessing the importance of numerous potentially relevant temperature metrics, we identified degree days as a valuable predictor of fish habitat, and predicted declines in walleye recruitment in spite of increased availability of thermal habitat when quantified using other metrics. Temperature influences fishes through multiple pathways (Magnuson et al., 1979) and the thermal niche of a species includes multiple life stages with differing thermal requirements (Hokanson, 1977; Wismer & Christie, 1987). A cumulative metric of temperature such as degree days may encompass myriad temperature influences in a single value and provide a measure of the metabolically relevant temperature experience of a fish that is difficult to quantify based on temperature optima. Furthermore, because random forest modeling allows for interactions among predictor variables that are not identified a priori, we identify lake characteristics that interact with the effect of temperature. These interactions allow us to identify refugia where coolwater species are likely to persist in spite of increased temperatures.

Importantly, our fish models are purely correlational and cannot identify the mechanism through which temperature influences walleye recruitment and the relative abundance of largemouth bass. We did not predict that most Wisconsin lakes will become too warm for walleye or other coolwater species to survive; in fact, by some metrics, thermal habitat for walleye was expected to increase. Still, temperature as measured by degree days was negatively correlated with walleye recruitment success, and positively correlated with largemouth bass abundance. Our results do not rule out the possibility of negative interactions between the two species; the true mechanism may be competition/predation between largemouth bass and walleye, or between each species and another, unmeasured variable that is itself correlated with degree days. We did not account for biotic interactions in our models, and such interactions can be important in determining species distributions (e.g., Guisan & Thuiller, 2005; Hein et al., 2012; Leathwick & Austin, 2001). However, abiotic variables are more important on a regional scale (Alofs & Jackson, 2015).

Numerous sources of uncertainty influence our results, and there are important limitations to our approach. We used a single emissions scenario for future projections (A2) that falls on the high end of all warming scenarios (Nakicenovic et al., 2000). As is the case for all studies attempting to predict future climates, we are uncertain about how much warming will actually occur. We believe using a high‐end warming scenario provides a useful benchmark for climate adaptation, because management plans that are robust to high‐end warming will likely be robust to less extreme futures. Furthermore, several studies using multiple scenarios showed little variability among scenarios compared to variability introduced by other modeling choices (e.g., GCMs or thresholds separating species presence from absence; Buisson et al., 2008; Hein et al., 2011; Thuiller, 2004). The three GCMs used here span most of the range of air temperature projections of all GCMs (Hostetler et al., 2011; Herb et al., 2014), and the variability among GCMs introduces further uncertainty regarding the rate, magnitude, and spatial variability of future climate change. An additional source of uncertainty is how lake temperatures will respond to changes in climate. We use a mechanistic model of lake temperatures that constrains temperatures based on thermodynamic and hydrodynamic first principles, which increases our confidence in using the model to project lake responses to temperatures warmer than what has been observed previously. This temperature model incorporates lake‐specific properties (such as size, depth, and water clarity) that contribute to heterogeneity among sites even when exposed to similar climatic conditions. As a result, these temperature projections are likely more robust than other studies over large extents and coarse resolutions (Guisan & Thuiller, 2005; Dormann, 2007). Still, like all models, our temperature model is a simplification of reality. For example, the thermodynamic model used here collapses within‐lake spatial variability into a single vertical axis. In reality, horizontal variability in water temperature may provide refugia for fish in some lakes.

Finally, predicted fish responses to warmer temperatures are uncertain. We converted probabilities of walleye recruitment and largemouth bass abundance to binary responses, and the use of somewhat arbitrarily defined thresholds adds uncertainty to our results (Thuiller, 2004; Wilson et al., 2005). For example, projected probabilities of the two species in lakes differ by a very small amount where we predict future ‘gains’ in walleye recruitment and ‘losses’ of high largemouth bass abundance, but these probabilities cross the threshold separating successful from failed recruitment or low from high abundance. Furthermore, projected fish relative abundance in the mid‐21st century differed substantially among GCMs (Table 3). Interestingly, even though late‐century water temperature predictions varied among GCMs, predicted fish relative abundance was relatively consistent due to the threshold responses of each species to water temperature (Fig. 1). After exceeding this threshold, any further increases in temperature were predicted to have minimal effects. However, these future projections step outside the realm of past observation in terms of water temperature degree days.

When conditions exceed the range of observed data, projections of species distributions based on correlations are uncertain, because relationships could change as climate changes (Dormann, 2007; Elith & Leathwick, 2009; Fitzpatrick & Hargrove, 2009). Behavioral responses, rapid evolution, or food web shifts can lead to surprises in species responses to climate change that are not predictable based on correlative species distribution models such as those employed here (e.g., Pearson & Dawson, 2003; Dormann et al., 2012). Alternative approaches are possible, where individual physiological responses to temperature changes are modeled mechanistically and scaled up to populations (e.g., Buckley et al., 2010); however, this approach is difficult to achieve on a scale and resolution that can influence management decision making (Guisan & Thuiller, 2005). Furthermore, predictions based on physiological responses to temperature can also be wrong if species exhibit unexpected behavioral shifts (e.g., Lawson et al., 2015). Although challenges exist regarding the application of correlative species distribution models such as those described here to future climate change scenarios, such models currently represent one of the only tools available for forecasting the impact of climate change on a wide range of species on management‐relevant scales (Guisan & Thuiller, 2005).

Ultimately, identifying the mechanisms behind observed correlations between species of interest and temperature may not matter for management decisions (Hilborn, 2016). Predicting lakes or lake types that are most likely to support desired species now and into the future will enable more effective resource allocation. Furthermore, by predicting the success of two recreationally important fish species simultaneously, managers can perhaps shift resources and target outreach efforts to focus on the species most likely to succeed in a given lake (in most cases, largemouth bass). We provide projected future responses to climate change of all managed lakes in Wisconsin in terms of walleye recruitment and largemouth bass relative abundance, and provide both lake‐specific probabilities and binary classes based on probability thresholds. Both are potentially useful for management. For lake‐specific management decisions and outreach with stakeholders, quantitative probabilities of the likelihood of supporting each species now and into the future may be most relevant. In contrast, big picture projections (e.g., Fig. 3) are likely to be most helpful for larger scale decisions such as allocation of resources for walleye and largemouth bass management on the statewide scale. For example, we identify dozens of lakes that are likely to support walleye recruitment even in a scenario of substantial warming, and management strategies designed to sustain and protect walleye could focus on these resilient lakes. Alternatively, managers can identify lakes where management actions such as walleye stocking will likely be the only option for providing a walleye fishery.

Understanding contemporary effects of climate and predicting the effects of future climate change on economically and ecologically important species is a critical need for resource management (Arvai et al., 2006; Heller & Zavaleta, 2009). To be most useful for management, such predictions should occur on the spatial scale and resolution at which decisions are made. Here, we used a process‐based model to predict lake‐specific responses to future climate conditions. In doing so, we accounted for heterogeneous responses of lakes to identical climatic conditions and avoided potential nonlinearities in the relationship between air temperature and water temperature (e.g., Mohseni & Stefan, 1999; Mohseni et al., 2002; Morrill et al., 2005) that could lead to overly pessimistic temperature projections. We also identified a single temperature metric (degree days) as the most influential of a large number of ecologically relevant temperature variables on two important fish species, suggesting that cumulative metrics of water temperature are more relevant to fishes than other metrics of ‘optimal’ thermal habitat (e.g., Wismer & Christie, 1987). Furthermore, we identify refuge lakes which are unlikely to be identified using models based solely on air temperature or those making predictions on coarse spatial resolutions. These local refugia are important candidates for protection and for future research on the characteristics that make them resilient to regional temperature increases. We argue that models accounting for such small‐scale heterogeneity in response to climate are critical for both understanding species responses to climate change and designing adaptation strategies for natural resource agencies across multiple scales (Hulme, 2005; Keppel & Wardell‐Johnson, 2012; Keppel et al., 2015).