Climate change drives widespread shifts in lake thermal habitat

Lake surfaces are warming worldwide, raising concerns about lake organism responses to thermal habitat changes. Species may cope with temperature increases by shifting their seasonality or their depth to track suitable thermal habitats, but these responses may be constrained by ecological interactions, life histories or limiting resources. Here we use 32 million temperature measurements from 139 lakes to quantify thermal habitat change (percentage of non-overlap) and assess how this change is exacerbated by potential habitat constraints. Long-term temperature change resulted in an average 6.2% non-overlap between thermal habitats in baseline (1978–1995) and recent (1996–2013) time periods, with non-overlap increasing to 19.4% on average when habitats were restricted by season and depth. Tropical lakes exhibited substantially higher thermal non-overlap compared with lakes at other latitudes. Lakes with high thermal habitat change coincided with those having numerous endemic species, suggesting that conservation actions should consider thermal habitat change to preserve lake biodiversity. Using measurements from 139 global lakes, the authors demonstrate how long-term thermal habitat change in lakes is exacerbated by species’ seasonal and depth-related constraints. They further reveal higher change in tropical lakes, and those with high biodiversity and endemism.

G lobal warming increases lake surface temperatures worldwide 1,2 , which strongly influences lake functioning, thermal structures and ecosystem processes 3-5 . As lakes warm, the available thermal habitat over specific temperature ranges can shrink or expand, with consequences for organisms, depending on their thermal tolerances 6-8 . In some cases, suitable thermal habitats may shrink or expand to the extent that native species become poorly adapted and non-native species thrive 9-11 . Changes in thermal habitat may be especially impactful in lakes because many species are ectothermic and, as on islands 12 and mountaintops 13 , they are partially restricted by the boundaries of lakes. Thus, the shrinkage and expansion of thermal habitats raise important concerns about how climate change affects lake ecosystems 14 and the implications of these changes for the threatened biodiversity that lake ecosystems currently support 15 .
With climate change, lakes are generally assumed to gain warm and lose cold thermal habitats. However, thermal habitat change in lakes is complex, as temperatures and temperature trends can vary vertically 16-18 , horizontally 19-21 and seasonally 21,22 within lakes. Furthermore, a substantial proportion of lakes exhibit cooling in their deeper waters, at least during stratified periods 16-18,23-25 , resulting in volumetrically cooler lakes in some regions 23 . This bottom-water cooling can be itself a thermal response of lakes to surface warming via a strengthening of thermal stratification 17,26 , which shields bottom waters from downward transmission of surface heating 25 . Furthermore, reductions in water clarity, such as temperature-induced increases in phytoplankton biomass 27 , can also cause deep-water cooling in lakes 28 . These contrasting mechanisms of surface and deep-water temperature change make it difficult to predict how thermal habitat may shift in response to observed global warming. Such knowledge is essential to improve our understanding of the vulnerability of lake ecosystems to climate change.
Some aquatic species may cope with a changing climate by shifting their seasonality (that is, phenology) or their depth distributions to track their suitable thermal habitat 6,29,30 . For example, some fishes with broad environmental tolerances may overcome local thermal habitat change by changing their depth and seasonality to take advantage of new resources and species interactions 31-34 . However, specialists whose seasonality and depth are constrained by species interactions, life history or resources such as light, nutrients and oxygen, may be more susceptible to thermal habitat change 35 . For example, Planktothrix rubescens, a filamentous cyanobacterium, is adapted to light and thermal stratification conditions which occur only at specific depths and seasons [36][37][38] . Similarly, some species of Daphnia, a common herbivorous zooplankton genus, partially rely on photoperiod as a cue for diapausing eggs to develop in spring, thereby limiting their capacity to track earlier phytoplankton blooms as lake temperatures change 39,40 . Anoxic zones may hinder . c-f, White dashed lines delineate expected habitats for three model taxa which are ecologically restricted in their depth (species 1 (d); a hypothetical low-light specialist phytoplankton), seasonality (species 2 (e); a spring migratory fish) and both depth and seasonality (species 3 (f); a diapausing benthic invertebrate). Thermal habitat is summarized by the volume occupying 0.1 °C wide temperature bins for the lake as a whole (c) and for each of the three model taxa (d-f); x and y axes in c-f differ to help visualize the extent of overlap between the past and present temperature distributions.
the ability of aerobic species to move deeper in lakes to track their suitable thermal habitat. Assessing whether habitat constraints exacerbate the risk of thermal habitat change would improve our understanding of the capacity of lake ecosystems to cope with climate variation.
Here we quantify the observed long-term thermal change in 139 lakes distributed across six continents-representing 69.3% of the Earth's surface freshwater habitat by water volume (70.6 × 10 3 km 3 out of the 101.8 × 10 3 km 3 global total) 41 . Using decades of lake temperature depth profiles, including more than 32 million total temperature observations, we calculated thermal non-overlap as the core metric of change. We defined thermal habitat change as the difference between recent lake temperatures (second half of each lake's time series) compared with an earlier baseline period (first half of each lake's time series). Thermal habitat change was quantified as the non-overlapped area of the two temperature distributions (recent and baseline) as a percentage of the combined area of those distributions, following an established method 42 . Temperature distributions were volume-weighted to best capture the volumetric habitat available for species. The resulting values of thermal non-overlap are a measure of relative thermal change standardized against temperature variation in the baseline period. To evaluate whether depth and seasonal habitat restrictions exacerbate thermal habitat changes, we recalculated thermal non-overlap over various potential restricted habitat ranges for each lake. These potential habitat restrictions provide a proxy for the broad range of possible species' capacities to cope with temperature change (Fig. 1). We used boosted regression trees (BRT) to model lake-to-lake variability Vänern, Sweden (a,c); Kinneret, Israel (b,d); Giles, United States (e,g); and Trout, United States (f,h) by the change in volume of lake water within 0.1 °C temperature bins as a percentage of each lake's total volume (c,d,g,h). The purple to yellow colour scale reinforces the temperature gradient on the x axis and is consistent across panels to facilitate the comparison across lakes. Lakes shown here were chosen for representativeness; the same plots for all 139 lakes are available at https://doi.org/10.5281/zenodo.4686874.
in thermal non-overlap to determine which lake characteristics are most strongly related to variability in thermal habitat change.
Over the past several decades, volume-weighted, whole-lake temperatures have increased in 77% of lakes with a mean trend of +0.12 °C per decade. However, whole-lake warming rates failed to reveal the complexity of thermal habitat change across each lake's temperature spectrum-most lakes exhibited serial losses and gains along their temperature gradients ( Fig. 2) with implications for which species would be most affected. After controlling for differences in the length and seasonal coverage of each lake's temperature time series and standardizing against a null estimate (Methods), we found that the mean thermal non-overlap across lakes from 1978 to 2013 was 6.2% and the distribution was skewed (median = 5.2%, inner quartile range = 4.5-6.9%; Fig. 3). Thus, 6.2% of the cumulative temperature distributions across both time periods are composed of either thermal habitat losses or gains over specific temperature ranges.
To evaluate whether habitat restrictions exacerbate the risk of thermal habitat change, we recalculated thermal non-overlap over various potential restricted habitat ranges for each lake. We considered two dimensions of habitat restrictions-depth and seasonality. The intensity of habitat range restrictions was scaled on a continuous gradient from 0 to 0.95 to facilitate comparisons across lakes, where 0 corresponds to no habitat restriction and 0.95 corresponds to habitats that are restricted to 5% of the available depths or days of the year. Thermal non-overlap was highest for the most restricted habitats (habitat restrictions of 0.95 in Fig. 4). When hypothetical habitats were restricted to 5% of all available depths, average thermal non-overlap increased from 6.2% to 9.7% (difference of 3.5%; Fig. 4a). When habitats were restricted to 5% of all available days of the year, average thermal non-overlap increased from 6.2% to 11.0% (difference of 4.8%; Fig. 4a). When habitats were restricted by both depth and season, strong synergistic interactions resulted in an increase in thermal non-overlap from 6.2% to 19.4% (difference of 13.2%). This interactive effect exceeded the additive effect when temperatures were compared at either restricted depths or restricted days of the year alone (13.2% > 3.5% + 4.8%).
The relative importance of seasonal versus depth habitat restrictions for thermal non-overlap depended on lake characteristicsseasonal habitat restrictions affected thermal non-overlap most strongly in shallower lakes (for example, Pesiöjärvi, Müggelsee and Annie; Fig. 4b-d) whereas depth habitat restrictions affected thermal non-overlap most strongly in deeper lakes (for example, Tahoe, Zürich and Tanganyika; Fig. 4e-g). Although we expect the most severe habitat restrictions (values of 0.95) to be rare in lakes, any organism with a more limited habitat range due to various ecological constraints would be likely to face greater thermal non-overlap. Conversely, generalist species whose seasonality or depth is unencumbered by strict ecological constraints would be markedly less likely to face thermal habitat changes.
The BRT analysis showed that patterns in the magnitude of thermal non-overlap among lakes were strongly associated with mean lake depth. Thermal non-overlap was higher in deeper lakes, especially when seasonal habitat restrictions were low (seasonal habitat restriction <0.65; Fig. 5). More intense seasonal habitat restrictions elevated thermal non-overlap values in shallow lakes, thereby minimizing the difference between thermal non-overlap in deep and shallow lakes (Fig. 5). This interaction may have been caused by the effect of baseline temperature variation on values of thermal non-overlap (Extended Data Fig. 1). Shallow lakes typically have larger seasonal temperature variation, which would tend to reduce thermal non-overlap. Thus, when seasonal habitat restrictions are applied to shallow lakes, the thermal non-overlap-enhancing effects of high seasonal temperature variation are minimized. In the light of this feature, we would expect the costs of various habitat restrictions in terms of thermal non-overlap to depend strongly on lake depth.
Patterns across lakes in the magnitude of thermal non-overlap were strongly associated with latitude (relative importance of all predictors shown in Extended Data Fig. 3). We expected temperate and arctic lakes to have greater thermal non-overlap because lake surface warming rates tend to be greater there 1 . However, in contrast to our expectation, tropical lakes exhibited substantially higher thermal non-overlap compared with lakes at other latitudes (Figs. 4 and 5). High thermal non-overlap in the tropics may arise from the thermal habitat change metric we used, which is affected by the baseline temperature variability (Extended Data Fig. 1). Smaller interand intra-annual temperature variation in tropical lakes upweights their magnitude of temperature change and produces relatively high thermal non-overlap. Quantifying thermal habitat change in this way may be more ecologically meaningful because the effects of temperature change on species depend on the breadth of species' thermal tolerances, which tend to be narrower when environmental temperature variation is low 43-45 . Overall, these analyses suggest that lake ecosystem sensitivity to climate change, as characterized by thermal non-overlap, differs considerably from the global pattern in the magnitude of lake surface warming. Despite slower surface warming rates, habitat changes are likely to be felt most strongly in tropical lakes where biodiversity may be most affected. The global pattern in thermal non-overlap may have differed from that of whole-lake warming rates (Extended Data Fig. 4) because transparency changes sometimes shield deep waters from downward transmission of heat 25 . Shielding would counterbalance calculations of whole-lake warming rates when integrating across the entire water column because some depths warm whereas others cool. When calculating whole-lake warming rates, lakes which exhibit deep-water cooling often have moderate or no significant differences in whole-lake mean annual temperature 25 . However, in our thermal non-overlap calculations, simultaneous warming and cooling over different parts of the water column would both elevate values of thermal non-overlap because the metric is agnostic toward the direction of change. Therefore, thermal non-overlap calculations reflect more than just heat content and, importantly, whole-lake temperature trends do not adequately describe the magnitude of thermal habitat change.
Some species may benefit from the changes in thermal habitat observed here. For instance, in Lake Zürich, the net growth rate of P. rubescens is maximized in warmer months (July-September) at specific depths (5-20 m) due to light, temperature and density stratification conditions. Restricting volumetric thermal habitat over this range increases the non-overlap value from 16% for the lake overall to 22% for the most suitable habitat for P. rubescens (Fig. 6). In response, P. rubescens growth is anticipated to increase in Lake Zürich 36,37 because warmer temperatures and higher thermal stratification are typically more suitable for the species' growth 46 . While some individual species with highly restricted habitats may benefit from thermal non-overlap 46 , thermal non-overlap is generally expected to increase the overall likelihood of species extinctions and community disruptions 14 . It has even been warned that climate change may drive certain species to local extirpations via changes to lake thermal habitat 47 . Our analysis provides a more global context for these concerns and goes further by showing that, in many cases, species will require the absence of habitat restrictions to effectively track their suitable thermal conditions. Shifts across depth and season are already widely observed in lakes for a variety of taxa 32,48,49 . However, the inability of most species to adopt these coping strategies 34,39,48 may increase the likelihood of community disruptions due to the disappearance of habitat over their suitable thermal ranges. Importantly, lakes with high thermal non-overlap closely coincide with those identified as containing a global heritage of freshwater biological diversity and endemism 50 , including Lakes Baikal (thermal non-overlap without habitat restrictions = 5.5%), Biwa (9.9%), Tanganyika (15.4%) and Victoria (18.5%). In these lakes, elevated risks of extinction due to thermal habitat change are possible, as is the disruption and disaggregation of lake ecosystems due to mismatched shifts across species 6,7,39 . Given the potential for the formation of new communities as thermal habitats change, the link between overall thermal habitat changes and ecological responses could fundamentally alter the nutrient and energy pathways that drive these ecosystems 6,7,39 .
Changes in lake thermal habitat may be especially impactful because most lake taxa are ectothermic, and many are partially restricted by the boundaries of lakes. Given these vulnerabilities to thermal habitat changes, standard conservation measures such as in-lake protected areas 51 may be ineffective at fully stemming the negative consequences of thermal habitat change 47 . Instead, lake conservation efforts focused on enhancing the capacity for organisms to shift across depth and season (for example, by reducing the size of lake anoxic zones) may help some species avoid thermal non-overlap. Furthermore, species that are active year-round in lakes may be sensitive to the full annual cycle of temperature without a chance to avoid seasons with unsuitable thermal habitat. Thus, for species habitats that are restricted across season or depth, conservation actions focused on enhancing within-lake shifts will be insufficient. Instead, efforts focused on enhancing connectivity among lakes (for example, dam removal) may be more effective at ameliorating the negative consequences of thermal habitat change 52 .
Lake temperature change can indirectly influence the habitat available for lake species aside from the direct temperature effect. For instance, lake temperature change influences underwater light availability and dissolved oxygen concentrations 28,53-55 . Anoxic zones may increase in extent or duration as lakes warm, further restricting the depths and seasons that can be occupied by aerobic organisms. Light penetration may increase or decrease with changes in water transparency driven by climate change with added consequences for the habitat available for photosynthetic organisms. Future studies focused on specific taxa should include these other important determinants of the habitat space available to them. The modular nature of the non-overlap calculations used herein allows for such   Fig. 6 | Thermal habitat change for the restricted habitat of P. rubescens in Lake Zürich. a, Epifluorescence microscopy photo of P. rubescens, a filamentous cyanobacterium present in Lake Zürich, with its optimal habitats occurring from July to September at depths of 5-20 m. b,c, Thermal habitat change (percentage of non-overlap) in this restricted range increases from 16% for Lake Zürich overall (b) to 22% for the habitat of P. rubescens (c). Credit: photo in a, Thomas Posch.
work incorporating a wider variety of environmental variables that determine habitat suitability 14 .
The redistribution of life on Earth is a major ecological response to anthropogenic global warming. As Earth warms, lake thermal habitats may shrink, expand or shift to seasons or depths where ecological interactions, life histories or resources limit species' growth and reproduction. These thermal shifts will inevitably have consequences for the species that lakes currently support. Forecasts of how species will respond to thermal non-overlap will always be less certain because they are often extrapolated from present conditions. However, with the results presented here, we have a framework for developing testable hypotheses relating lake thermal habitat change to biotic change over time for a large number of observed ecological time series. Given the risk of lake biodiversity loss, there is considerable societal value in resolving the connections between our results and real ecological responses in lakes.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41558-021-01060-3. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.

Overview.
We used long-term time series of lake temperature profiles to determine the magnitude of thermal habitat change in 139 widely distributed lakes. Time series were interpolated across depth and season to generate data with consistent resolutions across lakes. To assess temperature change, we used a metric, 'thermal non-overlap' , based on the percentage of two kernel density estimations of lake temperature which are non-overlapping. We calculated the metric for a range of plausible seasonal and depth habitat restrictions for aquatic species in the face of climate change. We used BRT to explain variability across lakes in their thermal habitat non-overlap as a function of lake characteristics (mean depth and latitude), characteristics of the time series for each lake (starting day of the year, ending day of the year, starting year and ending year, average number of sampling dates per year, long-term trend in the number of sampling dates per year, long-term trend in the yearly seasonal range of sampling dates), the habitat restriction values (season and depth) and the location of the time series delineation for thermal non-overlap calculations (30th, 50th and/or 70th quantiles of the years included in each lake's time series).

Study sites.
We compiled long-term lake temperature data from 139 lakes across the globe. Temperature variations in many of these lakes have already been linked to climate change 1,2,19,20,57,58 , but temperature change in at least one lake may be partially due to background climate variation in addition to anthropogenic climate change (Atlantic Multidecadal Oscillation in Lake Annie) 59 . The lakes included in our analysis represent a wide range of surface area (0.02 to 68,800 km 2 ), maximum depth (2.3 to 1,642 m), latitude (60 °S to 69 °N) and elevation (−212 to 1,987 m above sea level) (see Supplementary Table 1 for more information).
Temperature data. In total, we used more than 32 million lake temperature measurements for our analyses. The number of observations per lake ranged from 368 (Lake Stensjon) to 7,636,767 (Lake Superior) with approximately 232,000 observations per lake on average. Temperature data from each lake came from in situ temperature profiles 60-64 for lakes smaller than 169 km 2 and from a combination of in situ temperature profiles and remotely sensed surface water temperatures for 21 larger lakes. Remote sensing data were used in recognition that temperature and warming rates can vary substantially across latitude and longitude for large lakes 19-21 .
The mean length of the temperature time series was 36 years with a range from 15 to 101 years. All lakes had temperature data which started in the year 2000 or earlier and ended in 2000 or later. Lakes had on average 29 temperature profiles per year (inner quartile range: 7-26). In situ temperature data were measured using a wide variety of temperature sensors. Data collection methods included regularly collected discrete temperature profiles, high-resolution thermistor chains and other commonly accepted tools for measuring aquatic temperature. The in situ data are publicly available through the environmental data initiative 60 .
Remotely sensed lake surface temperatures were measured using the Advanced Very High-Resolution Radiometer (AVHRR) and processed by the Group for High Resolution Sea Surface Temperature (GHRSST) project 65 . AVHRR data have been validated against buoy data from the North American Great Lakes and found to have a root mean squared error of 0.55 °C compared with in situ measurements 2 . AVHRR temperature data were included to capture horizontal variability in temperature and warming in 21 of the 139 lakes that would not be captured by temperature profiles from a single central location [19][20][21] . AVHRR data were pooled with in situ data for temperature interpolation.
Temperature interpolation. Temperature data were spatially and temporally interpolated for each lake. All temperature profile data were first linearly interpolated across depth because temperature variability with depth is highly constrained by lake physics and typically allows for robust interpolations. The largest data gap over which depth interpolation occurred was 0.1 × mean depth of each lake. Following interpolation across depth, data were interpolated across time using standard spline interpolation models with a Kalman filter 66 . The model output was used to fill data gaps to produce a continuous, daily time series over the day of the year range for which temperature profiles had been regularly measured. Some times of the year were excluded from specific lakes because they lacked regular measurements throughout the length of the long-term time series. Thus, the same starting and ending day of the year was used for each lake throughout its time series, and was often shorter than the full annual cycle (Supplementary Table 1). The largest gap in time over which interpolation occurred was 30 days and this included extrapolations for lakes with missing data at the beginning or end of seasonal coverage in a specific year. Years with longer gaps were omitted from the analysis and the length of the seasonal coverage was optimized to minimize the number of years that needed to be removed. For large lakes with many sampling points (for example, Baikal, Superior, Victoria), temperature data were divided into 1,000 km 2 latitude-longitude bins and interpolated across depth and across time separately for each bin. The mean seasonal coverage of the interpolated lake time series was 245 days per year with a minimum of 17 days per year and a maximum of 365 days per year.
The interpolated temperature output had a daily temporal resolution and a depth resolution which varied continuously over depth. At the lake surface, we interpolated temperatures every 0.1 m (for example, 0 m, 0.1 m, 0.2 m), to every 1 m starting at a depth of 10 m (for example, 10 m, 11 m, 12 m) and every 100 m starting at a depth of 1,000 m (for example, 1,000 m, 1,100 m, 1,200 m). These depth increments were used because they consistently gave good coverage over all major lake strata, regardless of each lake's morphometric characteristics, while minimizing computational intensity by eliminating redundancy within lake strata.
Thermal habitat non-overlap calculations. After interpolating the temperature data across depth and season for each lake, we bisected it into an early part (part a) and a later part (part b). Parts a and b were iteratively delineated at three points positioned serially along the time series-at the 30th, 50th and 70th quantiles. We averaged the final non-overlap values across these three delineations for each lake so that the results depended less on the somewhat arbitrary decision of where to split the time series. For each delineation, we randomly sampled 10,000 temperature values from each of parts a and b. This was repeated ten times resulting in a total of 300,000 temperature values across all three time series delineations and all ten repetitions for each lake (10,000 × 3 × 10). The sampling probability for temperature values in each comparison was weighted by the volume increment associated with each temperature value (depth increment (I d ) × cross-sectional area at each depth (C d )). I d was calculated as the difference between the depth of the sampled temperature value and the next depth in the depth resolution of the interpolated temperatures. C d at each depth for each lake was calculated using standard, three-parameter models for estimating lake cross-sectional area based on surface area, maximum depth and mean depth 67 . For large lakes with temperature data at multiple locations across latitude and longitude, C d was divided by the number of latitude-longitude bins used for each lake. Temperature values from large lakes were sampled regardless of their associated latitude-longitude bins. As a result of the volume-weighting procedure, temperature measurements were sampled in proportion to the volume of water represented by each value, with temperatures representing larger volumes being sampled more often. As a consequence of this volume-weighting procedure, the resulting temperature distributions were robust to moderate changes in the depths used for the temperature interpolation ( Supplementary Fig. 1).
We defined thermal non-overlap (TNO) as the symmetric difference (Ө) between the kernel density estimations of temperature values from parts a and b of the time series as a proportion of the union (∪) of both kernel density estimations, following an established method 42 . Conversely, we defined the thermal habitat overlap (as opposed to non-overlap) as the intersection (∩) of the kernel density estimations as a proportion of the union (∪) of both distributions. All values were converted to percentages by multiplying by 100.
We used simulations to test the sensitivity of TNO to changes in mean and s.d. of temperature. We primed these simulations with three baseline temperature distributions all with a mean of 15 °C but with varying s.d. (4, 6, 8 °C). We simulated a range of additional temperature distributions by increasing and decreasing the mean and s.d. of the baseline temperature distributions and then calculated the corresponding values of TNO. The simulated change in both mean and s.d. varied from −3 to +3 °C. We found that TNO was sensitive to changes in mean and s.d. but was slightly more sensitive to reductions in s.d. compared with increases. TNO values also depended on the baseline s.d., such that lower starting s.d. elevates values of non-overlap given an equivalent change in temperature (Extended Data Fig. 1).
We also quantified null values of thermal non-overlap (TNO o ) by repeating the thermal non-overlap calculations but where parts a and b were defined by randomly dividing the individual years of data into two separate groups as opposed to sequentially dividing them along the time series.
To calculate standardized thermal non-overlap (TNO s ), we subtracted TNO o from TNO thereby setting the null expectation to zero.
In this case, if the temperature distributions in the recent and baseline time periods were identical, the TNO s would equal approximately zero. Values different from zero reflect a combination of random noise and long-term temperature change. All non-overlap values described in the main text and shown in Figs. 2-6 reflect values of TNO s . A comparison between raw values of TNO and TNO o can be found in Extended Data Fig. 5. Thermal non-overlap values and the null values were calculated using the 'overlap' function from the 'overlapping' package 42 in the R environment for statistical computing and visualization. In the function, we set the number of equally spaced points at which the overlapping kernel density estimation is evaluated to 100 for all comparisons because it minimized the values of TNO o (we considered a range of values from 5 to 10,000).
To assess the effect of seasonal habitat restrictions (S limit ) and volumetric habitat restrictions (V limit ), we modified equations (1)-(3) by comparing temperature values only from a specified range of depths and/or days of the year. We considered a range of habitat restrictions scaled from 0 to 0.95, where 0.95 is the most restrictive (temperature values were compared from within bins equivalent to 1/20th of the available seasonal and volumetric habitat) and 0 is the least restrictive (temperature values were compared regardless of season and depth). We focused our interpretations on the unitless habitat restrictions (scaled from 0 to 0.95) instead of in units of days or m 3 so that habitat restrictions could be more readily compared across lakes. Comparing a V limit value of 0.8 across lakes of different sizes assumes that a habitat restriction of 2 m 3 in a 10 m 3 lake would be comparable to a 20 m 3 habitat delineation in a 200 m 3 lake. The actual size of the seasonal habitat restrictions for each lake in units of days were calculated using the value of S limit as follows: where S is the seasonal habitat restriction in units of days, doy max is the maximum day of the year of the lakes' seasonal coverage, doy min is the minimum day of the year of the lakes' seasonal coverage and S limit is the seasonal habitat restriction scaled from 0 to 0.95. For example, in a lake with a seasonal coverage from day of the year 1 to day of the year 365, with an S limit value of 0.75, we compared randomly selected temperatures from time periods a and b separately for four seasonal bins (days of the year 1-91, 92-183, 184-273 and 274-365). Similarly, the actual size of the volumetric habitat restrictions (V) for each lake in units of m 3 were calculated using the value of V limit as follows: where V is the volumetric habitat restriction in units of m 3 , volume is the lake's total volume and V limit is the volumetric habitat restriction value scaled from 0 to 0.95. For example, if a lake with a volume of 100 m 3 had a V limit value of 0.75, we randomly selected temperature values from time periods a and b which were within four 25 m 3 (100 m 3 × (1 − 0.8)) bins. Volume bins were subsequently translated into sequential depth bins for the purpose of temperature value selection, making them functionally depth limits, and they are presented as such in the main text. We factorially combined a discrete series of values for S limit and V limit (0, 1/2, 2/3, 5/6, 8/9, 12/13 and 19/20) to test a range of combined seasonal and volumetric habitat restrictions that do not require the overlap or truncation of bins. For reference, habitat restrictions are presented visually for hypothetical 'Species 1' (S limit = 0, V limit = 0.8), 'Species 2' (S limit = 0.8, V limit = 0) and 'Species 3' (S limit = 0.8, V limit = 0.8) examples (Fig. 1). These limits reflect hypothetical restrictions in a species' habitat due to ecological factors and approximate the habitat available for a low-light specialist phytoplankton (species 1), a spring migratory fish (species 2) and a diapausing benthic invertebrate (species 3). In Fig. 6, the species habitat restriction values for P. rubescens were S limit = 0.74, V limit = 0.89 (Fig. 6).
Explaining variability in thermal habitat non-overlap. We used BRT to explain lake-to-lake variability in thermal habitat change (percentage of non-overlap) while accounting for differences in the temporal coverage of each lake's time series. The predictor variables in the BRT were the starting year of the time series, ending year of the time series, starting day of the year of the seasonal coverage, ending day of the year of the seasonal coverage, average number of sampling dates per year, linear trend (Theil-Sen slope) in the average number of sampling dates per year, linear trend (Theil-Sen slope) in the yearly extent of the time series' seasonal coverage, lake mean depth, absolute latitude (degrees from the Equator), seasonal habitat restriction, depth habitat restriction and time series delineation. Geospatial and morphometric data for each lake is available from the previously published HydroLAKES database 41 . Of the available lake characteristics, we used latitude and mean depth because they were most strongly correlated to TNO s values and because they were least correlated to the other predictors in the model. We used a 100-fold cross-validation with a 70-30% split by lake (that is, 70% of lakes were used in each BRT). Model results were averaged to ensure that the patterns described therein were robust to the exclusion of some lakes. We optimized the learning rate for each BRT by iteratively running the model with smaller and smaller learning rates (from 0.8, 0.4, 0.2, 0.1, 0.05 to 0.025) until the number of trees in the model was greater than 1,000, as suggested in previous literature 68 . We found that the BRT performed well in cross-validation-the correlation between predicted and observed values in the test datasets from the 100-fold cross-validation was moderate on average across models (r = 0.56, Kendall's rank correlation; see full goodness-of-fit summary statistics in Extended Data Fig. 6). The correlation between the predicted and the observed values was high (r = 0.76, Kendall's rank correlation) when predictions were averaged across BRT. We found minimal patterning in the model residuals when comparing the model residuals with each predictor variable used in the BRT (Extended Data Fig. 7).
To calculate lake-specific mean thermal non-overlap values and facilitate comparison across lakes, we used the BRT to remove the variation in thermal non-overlap attributable to the starting year of the time series, ending year of the time series, starting day of the year of the seasonal coverage, ending day of the year of the seasonal coverage, average number of sampling dates per year, linear trend (Theil-Sen slope) in the average number of sampling dates per year and the linear trend (Theil-Sen slope) in the yearly extent of the time series' seasonal coverage of each lake's time series, following previously published work 24 . We did this by setting the values for these variables to their median and using the BRT to make a prediction for each lake with these medians as predictors, along with each lake's observed values for mean depth, absolute latitude, seasonal habitat restriction, depth habitat restriction and time series delineation. The residuals from the BRT were then added back to the predicted values used in further analyses and plotting. The mean lake-specific thermal dissimilarities were calculated as the average across all seasonal habitat restrictions (S limit ), depth habitat restrictions (V limit ) (0, 1/2, 2/3, 5/6, 8/9, 12/13 and 19/20) and all three time series delineations. The statistical significance of these lake-specific thermal non-overlap values was estimated on a continuous gradient and calculated using a Wilcoxon signed-rank test. In the test, we compared TNO values to TNO o values separately for each combination of time series delineation, seasonal habitat restriction and depth habitat restriction (n = 108). The average P values from these tests for each lake are shown in Supplementary Table 1.
We compared thermal non-overlap values to a more widely used metric of whole-lake thermal change-whole-lake temperature trends. Whole-lake temperature trends were calculated based on the annual averages of all temperature values sampled for the pairwise thermal non-overlap calculations to maximize the comparability of the resulting temperature trends and thermal non-overlap values. Due to the temperature sampling probability being volume-weighted, the temperature trend was also indirectly volume-weighted. Temperature trends were calculated using Theil-Sen slopes applied to annual mean temperatures and the statistical significance of each trend (P value) was calculated using a bootstrapped one sample Wilcoxon signed-rank test with 1,000 repetitions. The input data for the Wilcoxon signed-rank test were the complete list of all slopes derived from all pairwise combinations of points in the time series. The number of pairwise slopes used in each repetition of the Wilcoxon signed-rank test was equal to the number of years of temperature data for each lake. Whole-lake temperature trends and thermal non-overlap values were not strongly correlated (r = 0.10, Kendall's rank correlation coefficient; Extended Data Fig. 4). All statistics and graphics were produced in the R statistical computing environment 69 . Fig. 1 | Sensitivity of thermal habitat change values to changes in mean and standard deviation of three hypothetical temperature distributions. The three hypothetical baseline temperature distributions considered here had a mean of 15 °C and standard deviations of 4, 6, and 8 °C (a). We simulated the effects of temperature changes (increases and decreases in both the mean and the standard deviation of the baseline temperature distributions) on resulting values of thermal habitat change (% non-overlap). The effects on thermal habitat change were simulated for each of the three baseline conditions (b-d). Values of thermal habitat change (% non-overlap) were sensitive to changes in mean and standard deviation of temperature as well as to the standard deviation of the baseline temperature distribution. Higher standard deviation of the baseline distribution lead to lower values of thermal habitat change. Fig. 2 | Partial dependency plots for each of the 12 predictors in the boosted regression trees. Colored lines are the raw partial dependence plot values from each of the 100 different BRTs which were fit to random subsets (70%) of the 139 lakes included in this analysis and the black lines are the locally weighted scatterplot smoothed (LOESS) lines across all BRTs. The panels are ordered by the relative importance of each predictor in the model (a-l) with the most important predictors in the upper left (a) and the least important in the bottom right (l). Note that the scales on the x-and y-axes vary from plot to plot.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Software and code
Policy information about availability of computer code Data collection No software was used to collect the lake temperature data.

Data analysis
The R statistical computing environment was used for the analysis of all data and for the production of all figures in this manuscript. All code can be found under the identifier, DOI: 10.5281/zenodo.3894908 For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability The remote sensing lake temperature data used here can be found under the identifier doi.org/10.5067/GHAAO-4BC02 and through the Earthdata website (https:// earthdata.nasa.gov/). The in situ lake temperature data used here are available through the Environmental Data Initiative (EDI) data portal (https:// portal.edirepository.org/nis/home.jsp#) under the identifier, doi:10.6073/pasta/f03d7d682eae2f467642e4260686ea15. Geospatial and morphometric data for each lake are available from the HydroLAKES database under the identifier, doi: 10.1038/ncomms13603 and can be found at http://www.hydrosheds.org.

nature research | reporting summary
April 2020 Field-specific reporting Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

Life sciences Behavioural & social sciences Ecological, evolutionary & environmental sciences
For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf

Ecological, evolutionary & environmental sciences study design
All studies must disclose on these points even when the disclosure is negative.

Study description
We use long-term time series of lake temperature profiles to determine the magnitude of thermal habitat change in 139 widely distributed lakes. Time series were interpolated across depth and across seasons to generate data with consistent resolutions across lakes. To assess temperature change, we used a metric, "thermal novelty", based on standardized Euclidean distances applied to the full length of each lake's interpolated time series. We calculated the metric for a range of plausible seasonal and depth shift capacities for aquatic species in the face of climate change. We used a boosted regression tree (BRT) to explain variability across lakes in their thermal habitat novelty using lake characteristics (shore development index, residence time, mean depth, surface area, latitude, longitude, and elevation) while controlling for the characteristics of the time series for each lake (starting day of the year, ending day of the year, starting year, and ending year) and the shift limit (seasonal shift limit and depth shift limit).

Research sample
In situ temperature measurements were recorded digitally or in writing at the time of measurement.

Sampling strategy
Long-term lake water temperature data were assembled from as many monitoring stations around the world as possible.

Data collection
The data reported here reflect the collective efforts of myriad dedicated field crews, laboratory staff, data management and quality control staff, analysts and many others from a wide variety of nations, states, tribes, agencies, universities, and other organizations.
Timing and spatial scale The timing and spatial scale of temperature measurements are complex and vary from lake to lake. The duration and seasonal coverage for each the temperature profile time series from each lake are summarized in the supplementary material.

Data exclusions
No data were excluded.

Reproducibility
There are no experimental findings in our study.

Randomization
The lake water temperature data reported here were collected on an opportunistic basis from as many lakes as possible. The length and seasonal coverage of lake water temperature time series varied from lake to lake and this variation was accounted for using boosted regression trees.