<figure><img src="image_1.png" /></figure>

Check for updates

OPEN

# Climate change drives widespread shifts in lake thermal habitat

Benjamin M. Kraemer ID 1✉, Rachel M. Pilla ID 2, R. lestyn Woolway ID 3,4, Orlane Anneville ID 5, Syuhei Ban ID 6, William Colom-Montero ID 7, Shawn P. Devlin ID 8, Martin T. Dokulil ID 9, Evelyn E. Gaiser ID 10, K. David Hambright ID 11, Dag O. Hessen ID 12, Scott N. Higgins ID 13, Klaus D. Jöhnk ID 14, Wendel Keller 15, Lesley B. Knoll ID 16, Peter R. Leavitt ID 17,18, Fabio Lepori ID 19, Martin S. Luger ID 20, Stephen C. Maberly ID 21, Dörthe C. Müller-Navarra ID 22, Andrew M. Paterson ID 23,24, Donald C. Pierson ID 25, David C. Richardson ID 26, Michela Rogora ID 27, James A. Rusak ID 23,24, Steven Sadro ID 28, Nico Salmaso ID 29, Martin Schmid ID 30, Eugene A. Silow ID 31, Ruben Sommaruga ID 32, Julio A. A. Stelzer ID 1,33,34, Dietmar Straile ID 35, Wim Thiery ID 36,37, Maxim A. Timofeyev ID 31, Piet Verburg ID 38, Gesa A. Weyhenmeyer ID 25 and Rita Adrian 1,34

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.

Global warming increases lake surface temperatures worldwide<fnref n="1" />,², which strongly influences lake functioning, thermal structures and ecosystem processes³⁻⁵. As lakes warm, the available thermal habitat over specific temperature ranges can shrink or expand, with consequences for organisms, depending on their thermal tolerances⁶⁻⁸. In some cases, suitable thermal habitats may shrink or expand to the extent that native species become poorly adapted and non-native species thrive⁹⁻¹¹. Changes in thermal habitat may be especially impactful in lakes because many species are ectothermic and, as on islands¹² and mountaintops¹³, 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¹⁴ and the implications of these changes for the threatened biodiversity that lake ecosystems currently support¹⁵.

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¹⁶⁻¹⁸, horizontally¹⁹⁻²¹ and seasonally²¹,²² within lakes. Furthermore, a substantial proportion of lakes exhibit cooling in their deeper waters, at least during stratified periods¹⁶⁻¹⁸,²³⁻²⁵, resulting in volumetrically cooler lakes in some regions²³. This bottom-water cooling can be itself a thermal response of lakes to surface warming via a strengthening of thermal stratification¹⁷,²⁶,

<footnote n="1">¹Leibniz Institute for Freshwater Ecology and Inland Fisheries, Berlin, Germany. ²Miami University, Oxford, OH, USA. ³Centre for Freshwater and Environmental Studies, Dundalk Institute of Technology, Dundalk, Ireland. ⁴European Space Agency Climate Office, ECSAT, Harwell Campus, Didcot, UK. ⁵National Research Institute for Agriculture, Food and Environment (INRAE), University Savoie Mont Blanc, CARRTEL, Thonon-Les-Bains, France. ⁶University of Shiga Prefecture, Hikone, Japan. ⁷Uppsala University, Norrtälje, Sweden. ⁸University of Montana, Polson, MT, USA. ⁹University of Innsbruck, Mondsee, Austria. ¹⁰Florida International University, Miami, FL, USA. ¹¹University of Oklahoma, Norman, OK, USA. ¹²University of Oslo, Oslo, Norway. ¹³ISD Experimental Lakes Area Inc, Winnipeg, Manitoba, Canada. ¹⁴Commonwealth Scientific and Industrial Research Organization, Canberra, Australian Capital Territory, Australia. ¹⁵Laurentian University, Sudbury, Ontario, Canada. ¹⁶University of Minnesota, Lake Itasca, MN, USA. ¹⁷University of Regina, Regina, Saskatchewan, Canada. ¹⁸Queen's University Belfast, Belfast, UK. ¹⁹University of Applied Sciences and Arts of Southern Switzerland, Canobbio, Switzerland. ²⁰Federal Agency for Water Management, Mondsee, Austria. ²¹UK Centre for Ecology & Hydrology, Lancaster, UK. ²²University of Hamburg, Hamburg, Germany. ²³Ontario Ministry of the Environment Conservation and Parks, Dorset, Ontario, Canada. ²⁴Queen's University, Kingston, Ontario, Canada. ²⁵Uppsala University, Uppsala, Sweden. ²⁶State University of New York at New Paltz, New Paltz, NY, USA. ²⁷CNR Water Research Institute (IRSA), Verbania Pallanza, Italy. ²⁸University of California-Davis, Davis, CA, USA. ²⁹Research and Innovation Centre, Fondazione Edmund Mach, San Michele all'Adige, Italy. ³⁰Eawag: Swiss Federal Institute of Aquatic Science and Technology, Kastanienbaum, Switzerland. ³¹Irkutsk State University, Irkutsk, Russia. ³²Department of Ecology, University of Innsbruck, Innsbruck, Austria. ³³University of Geneva, Geneva, Switzerland. ³⁴Freie Universität Berlin, Berlin, Germany. ³⁵University of Konstanz, Konstanz, Germany. ³⁶ETH Zürich, Zürich, Switzerland. ³⁷Vrije Universiteit Brussel, Brussels, Belgium. ³⁸National Institute of Water and Atmospheric Research, Hamilton, New Zealand. ✉e-mail: ben.m.kraemer@gmail.com</footnote>
a

<figure><img src="image_2.png" /></figure>

b

<figure><img src="image_3.png" /></figure>

c

<figure><img src="image_4.png" /></figure>

d

<figure><img src="image_5.png" /></figure>

e

<figure><img src="image_6.png" /></figure>

f

<figure><img src="image_7.png" /></figure>

**Fig. 1** | Simulated effects of lake warming on the exposure of hypothetical aquatic species to thermal habitat change. a,b, Simulated lake temperatures for a model lake with an average temperature increase of 1.14 °C between baseline (a) and recent (b) time intervals, consistent with known warming trends (1980s to 2010s)¹⁷,²⁴,⁵⁶. 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.

which shields bottom waters from downward transmission of surface heating²⁵. Furthermore, reductions in water clarity, such as temperature-induced increases in phytoplankton biomass²⁷, can also cause deep-water cooling in lakes²⁸. 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⁶,²⁹,³⁰. 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³¹⁻³⁴. 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³⁵. For example, *Planktothrix rubescens*, a filamentous cyanobacterium, is adapted to light and thermal stratification conditions which occur only at specific depths and seasons³⁶⁻³⁸. 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³⁹,⁴⁰. Anoxic zones may hinder
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

**Fig. 2** | Long-term changes in cumulative annual thermal habitat in lakes by volume. **a–h.** Volumetric thermal habitats are summarized for Lakes 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.

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³ km³ out of the 101.8 × 10³ km³ global total)⁴¹. 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⁴². 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
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

**Fig. 3** | Global patterns in lake thermal habitat change (percentage of non-overlap; 1978–2013). **a**, Map of thermal non-overlap where the variation attributable to the length and seasonal coverage of each lake's time series has been accounted for and removed from the values to facilitate comparisons across lakes. **b**, Thermal non-overlap values are shown relative to the cases where habitats are most restricted (habitats limited to 5% of the depth range and 5% of the seasonal range). The colour scales in **a** and **b** are assigned by each lake's rank in thermal non-overlap.

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 characteristics—seasonal 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
<figure><img src="image_18.png" /></figure>

**Fig. 4** | Thermal non-overlap as a function of the severity of habitat restrictions. a-g. The average thermal non-overlap across all lakes (a), and for individual lakes (b-g). Habitat restriction values are hypothetical constraints on any species' ability to move across depth or seasons to track their thermal habitat. Habitat restriction values are quantified in each plot as a proportion of the maximum from 0 to 0.95, where 0 is a habitat encompassing the entire depth range and seasonal extent and 0.95 is a habitat restricted to 5% of available depths (y axes) or days of the year (x axes). Lakes shown here were chosen for representativeness; the same plots for all 139 lakes are available at https://doi.org/10.5281/zenodo.4686874.

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¹. 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 inter- and 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
a

<figure><img src="image_19.png" /></figure>

b

<figure><img src="image_20.png" /></figure>

c

<figure><img src="image_21.png" /></figure>

d

<figure><img src="image_22.png" /></figure>

e

<figure><img src="image_23.png" /></figure>

f

<figure><img src="image_24.png" /></figure>

**Fig. 5** | Drivers of lake-to-lake variability in thermal non-overlap based on BRT analysis. **a,b,** Kernel density estimates for the lake mean depth (a) and lake latitude (b) across all 139 lakes. **c–f,** Coloured lines show thermal non-overlap values for different severities of seasonal and depth habitat restrictions as a proportion of the maximum from 0 to 0.95, where 0 is a habitat encompassing the entire depth range and seasonal extent of each lake and 0.95 is a habitat restricted to 5% of the available depths or days of the year. Lines represent 0.05-unit increments of habitat restriction severity. **a,b** share the x axes from **c–f**. The relationship between all 12 model predictors and the corresponding non-overlap values are shown in Extended Data Fig. 2.

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<sup>43-45</sup>. 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<sup>25</sup>. 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<sup>25</sup>. 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<sup>36,37</sup> because warmer temperatures and higher thermal stratification are typically more suitable for the species' growth<sup>46</sup>. While some individual species with highly restricted habitats may benefit from thermal non-overlap<sup>46</sup>, thermal non-overlap is generally expected to increase the overall likelihood of species extinctions and community disruptions<sup>14</sup>. It has even been warned that climate change may drive certain species to local extirpations via changes to lake thermal habitat<sup>47</sup>. Our analysis provides a more global context for these concerns and goes further by showing that, in many cases,
<figure><img src="image_25.png" /><figcaption>**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.</figcaption></figure>

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³²,⁴⁸,⁴⁹. However, the inability of most species to adopt these coping strategies³⁴,³⁹,⁴⁸ 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⁵⁰, 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⁶,⁷,³⁹. 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⁶,⁷,³⁹.

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⁵¹ may be ineffective at fully stemming the negative consequences of thermal habitat change⁴⁷. 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⁵².

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²⁸,⁵³⁻⁵⁵. 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
work incorporating a wider variety of environmental variables that determine habitat suitability¹⁴.

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.

Received: 25 May 2020; Accepted: 28 April 2021;
Published online: 3 June 2021

## References

1. O'Reilly, C. M. et al. Rapid and highly variable warming of lake surface waters around the globe. *Geophys. Res. Lett.* **42**, 10773-10781 (2015).

2. Schneider, P. & Hook, S. J. Space observations of inland water bodies show rapid surface warming since 1985. *Geophys. Res. Lett.* **37**, 22405-22409 (2010).

3. Comte, L. & Olden, J. D. Climatic vulnerability of the world's freshwater and marine fishes. *Nat. Clim. Change* **7**, 718-722 (2017).

4. Blois, J. L., Zarnetske, P. L., Fitzpatrick, M. C. & Finnegan, S. Climate change and the past, present, and future of biotic interactions. *Science* **341**, 499-504 (2013).

5. Woolway, R. I. et al. Global lake responses to climate change. *Nat. Rev. Earth Environ.* **1**, 388-403 (2020).

6. Tunney, T. D., McCann, K. S., Lester, N. P. & Shuter, B. J. Effects of differential habitat warming on complex communities. *Proc. Natl Acad. Sci. USA* **111**, 8077-8082 (2014).

7. McMeans, B. C. et al. The adaptive capacity of lake food webs: from individuals to ecosystems. *Ecol. Monogr.* **86**, 4-19 (2016).

8. Cline, T. J., Bennington, V. & Kitchell, J. F. Climate change expands the spatial extent and duration of preferred thermal habitat for Lake Superior fishes. *PLoS ONE* **8**, e62279 (2013).

9. Till, A., Rypel, A. L., Bray, A. & Fey, S. B. Fish die-offs are concurrent with thermal extremes in north temperate lakes. *Nat. Clim. Change* **9**, 637-641 (2019).

10. Jonsson, T. & Setzer, M. A freshwater predator hit twice by the effects of warming across trophic levels. *Nat. Commun.* **6**, 5992 (2015).

11. Muhlfeld, C. C. et al. Invasive hybridization in a threatened species is accelerated by climate change. *Nat. Clim. Change* **4**, 620-624 (2014).

12. Bellard, C., Russell, J., Hoffmann, B. D., Leclerc, C. & Courchamp, F. Adapting island conservation to climate change. Response to Andréfouët et al. *Trends Ecol. Evol.* **30**, 2-3 (2015).

13. Nogués-Bravo, D., Araújo, M. B., Errea, M. P. & Martínez-Rica, J. P. Exposure of global mountain systems to climate warming during the 21st century. *Glob. Environ. Change* **17**, 420-428 (2007).

14. Williams, J. W., Jackson, S. T. & Kutzbach, J. E. Projected distributions of novel and disappearing climates by 2100 AD. *Proc. Natl Acad. Sci. USA* **104**, 5738-5742 (2007).

15. Abell, R. et al. Freshwater ecoregions of the world: a new map of biogeographic units for freshwater biodiversity conservation. *BioScience* **58**, 403-414 (2008).

16. Winslow, L. A., Read, J. S., Hansen, G. J. A. & Hanson, P. C. Small lakes show muted climate change signal in deep-water temperatures. *Geophys. Res. Lett.* **42**, 355-361 (2014).

17. Kraemer, B. M. et al. Morphometry and average temperature affect lake stratification responses to climate change. *Geophys. Res. Lett.* **42**, 4981-4988 (2015).

18. Niedrist, G. H., Psenner, R. & Sommaruga, R. Climate warming increases vertical and seasonal water temperature differences and inter-annual variability in a mountain lake. *Climatic Change* **151**, 473-490 (2018).

19. Mason, L. A. et al. Fine-scale spatial variation in ice cover and surface temperature trends across the surface of the Laurentian Great Lakes. *Climatic Change* **138**, 71-83 (2016).

20. Woolway, R. I. & Merchant, C. J. Intralake heterogeneity of thermal responses to climate change: a study of large Northern Hemisphere lakes. *J. Geophys. Res. Atmos.* **6**, 3087-3098 (2018).

21. Kraemer, B. M. et al. Century-long warming trends in the upper water column of Lake Tanganyika. *PLoS ONE* **10**, e0132490 (2015).

22. Winslow, L. A., Read, J. S., Hansen, G. J. A. A., Rose, K. C. & Robertson, D. M. Seasonality of change: summer warming rates do not fully represent effects of climate change on lake temperatures. *Linnol. Oceanogr.* **62**, 2168-2178 (2017).

23. Tanentzap, A. J. et al. Cooling lakes while the world warms: effects of forest regrowth and increased dissolved organic matter on the thermal regime of a temperate, urban lake. *Linnol. Oceanogr.* **53**, 404-410 (2008).

24. Richardson, D. C. et al. Transparency, geomorphology and mixing regime explain variability in trends in lake temperature and stratification across Northeastern North America (1975-2014). *Water* **9**, 442 (2017).

25. Bartosiewicz, M. et al. Hot tops, cold bottoms: synergistic climate warming and shielding effects increase carbon burial in lakes. *Linnol. Oceanogr. Lett.* **4**, 132-144 (2019).

26. Woolway, R. I. & Merchant, C. J. Worldwide alteration of lake mixing regimes in response to climate change. *Nat. Geosci.* **12**, 271-276 (2019).

27. Ho, J. C., Michalak, A. M. & Pahlevan, N. Widespread global increase in intense lake phytoplankton blooms since the 1980s. *Nature* **574**, 667-670 (2019).

28. Rose, K., Winslow, L. & Read, J. Climate-induced warming of lakes can be either amplified or suppressed by trends in water clarity. *Linnol. Oceanogr. Lett.* **1**, 44-53 (2016).

29. Bartley, T. J. et al. Food web rewiring in a changing world. *Nat. Ecol. Evol.* **3**, 345-354 (2019).

30. Cheung, W., Close, C., Lam, V., Watson, R. & Pauly, D. Application of macroecological theory to predict effects of climate change on global fisheries potential. *Mar. Ecol.* **365**, 187-197 (2008).

31. Grieve, B. D., Curchitser, E. N. & Rykaczewski, R. R. Range expansion of the invasive lionfish in the Northwest Atlantic with climate change. *Mar. Ecol. Prog. Ser.* **546**, 225-237 (2016).

32. Hovel, R. A., Carlson, S. M. & Quinn, T. P. Climate change alters the reproductive phenology and investment of a lacustrine fish, the three-spine stickleback. *Glob. Chang. Biol.* **6**, 2308-2320, https://doi.org/10.1111/gcb.13531 (2017).

33. Campana, S. E. et al. Arctic freshwater fish productivity and colonization increase with climate warming. *Nat. Clim. Change* **10**, 428-433 (2020).

34. Dahlke, F. T., Wohlrab, S., Butzin, M. & Pörtner, H. O. Thermal bottlenecks in the life cycle define climate vulnerability of fish. *Science* **369**, 65-70 (2020).

35. Cheung, W. et al. Shrinking of fishes exacerbates impacts of global ocean changes on marine ecosystems. *Nat. Clim. Change* **3**, 254-258 (2012).

36. Gallina, N., Beniston, M. & Jacquet, S. Estimating future cyanobacterial occurrence and importance in lakes: a case study with *Planktothrix rubescens* in Lake Geneva. *Aquat. Sci.* **79**, 249-263 (2017).

37. Yankova, Y., Villiger, J., Pernthaler, J., Schanz, F. & Posch, T. Prolongation, deepening and warming of the metaliminion change habitat conditions of the harmful filamentous cyanobacterium *Planktothrix rubescens* in a prealpine lake. *Hydrobiologia* **776**, 125-138 (2016).

38. Dokulil, M. T. & Teubner, K. Deep living *Planktothrix rubescens* modulated by environmental constraints and climate forcing. *Hydrobiologia* **698**, 29-46 (2012).

39. Winder, M. & Schindler, D. E. Climate change uncouples trophic interactions in an aquatic ecosystem. *Ecology* **85**, 2100-2106 (2004).

40. Stross, R. G. Light and temperature requirements for diapause development and release in *Daphnia*. *Ecology* **47**, 368-374 (1966).

41. Messenger, M. L., Lehner, B., Grill, G., Nedeva, I. & Schmitt, O. Estimating the volume and age of water stored in global lakes using a geo-statistical approach. *Nat. Commun.* **7**, 13603 (2016).

42. Pastore, M. & Calcagni, A. Measuring distribution similarities between samples: a distribution-free overlapping index. *Front. Psychol.* **10**, 1089 (2019).

43. Tewksbury, J. J., Huey, R. B. & Deutsch, C. A. Ecology. Putting the heat on tropical animals. *Science* **320**, 1296-1297 (2008).

44. Deutsch, C. et al. Impacts of climate warming on terrestrial ectotherms across latitude. *Proc. Natl Acad. Sci. USA* **105**, 6668-6672 (2008).

45. Pinsky, M. L., Eikeset, A. M., McCauley, D. J., Payne, J. L. & Sunday, J. M. Greater vulnerability to warming of marine versus terrestrial ectotherms. *Nature* **569**, 108-111 (2019).
46. Posch, T., Köster, O., Salcher, M. M. & Pernthaler, J. Harmful filamentous cyanobacteria favoured by reduced water turnover with lake warming. *Nat. Clim. Change* 2, 809–813 (2012).

47. Heino, J., Virkkala, R. & Toivonen, H. Climate change and freshwater biodiversity: detected patterns, future trends and adaptations in northern regions. *Biol. Rev.* **84**, 39–54 (2009).

48. Adrian, R., Wilhelm, S. & Gerten, D. Life-history traits of lake plankton species may govern their phenological response to climate warming. *Glob. Chang. Biol.* **12**, 652–661 (2006).

49. Thackeray, S. J., Jones, I. D. & Maberly, S. C. Long-term change in the phenology of spring phytoplankton: species-specific responses to nutrient enrichment and climatic change. *J. Ecol.* **96**, 523–535 (2008).

50. Vadeboncoeur, Y., McIntyre, P. B. & Vander Zanden, M. J. Borders of biodiversity: life at the edge of the world's large lakes. *BioScience* **61**, 526–537 (2011).

51. McCullough, I. M., Skaff, N. K., Soranno, P. A. & Cheruvelil, K. S. No lake left behind: how well do U.S. protected areas meet lake conservation targets? *Linnol. Oceanogr. Lett.* **4**, 183–192 (2019).

52. O'Connor, J. E., Duda, J. J. & Grant, G. E. 1000 dams down and counting: dam removals are reconnecting rivers in the United States. *Science* **348**, 496–497 (2015).

53. Van Bocxlaer, B., Schultheiß, R., Plisnier, P.-D. & Albrecht, C. Does the decline of gastropods in deep water herald ecosystem change in Lakes Malawi and Tanganyika? *Freshw. Biol.* **57**, 1733–1744 (2012).

54. Palmer, M. E., Yan, N. D. & Somers, K. M. Climate Change drives coherent trends in physics and oxygen content in North American lakes. *Climatic Change* **124**, 285–299 (2014).

55. Kraemer, B. M., Mehner, T. & Adrian, R. Reconciling the opposing effects of warming on phytoplankton biomass in 188 large lakes. *Sci. Rep.* **7**, 10762 (2017).

56. Pilla, R. M. et al. Deeper waters are changing less consistently than surface waters in a global analysis of 102 lakes. *Sci. Rep.* **10**, 20514 (2020).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_26.png" /></figure>

**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/.

© The Author(s) 2021
## Methods

**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<sup>1,2,19,20,57,58</sup>, 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)<sup>59</sup>. The lakes included in our analysis represent a wide range of surface area (0.02 to 68,800 km<sup>2</sup>), 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<sup>60–64</sup> for lakes smaller than 169 km<sup>2</sup> 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<sup>19–21</sup>.

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<sup>60</sup>.

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 (GHRSSST) project<sup>65</sup>. 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<sup>2</sup>. 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<sup>19–21</sup>. 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<sup>66</sup>. 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<sup>2</sup> 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<sup>67</sup>. 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 ($\Theta$) between the kernel density estimations of temperature values from parts a and b of the time series as a proportion of the union (U) of both kernel density estimations, following an established method<sup>42</sup>. Conversely, we defined the thermal habitat overlap (as opposed to non-overlap) as the intersection ($\cap$) of the kernel density estimations as a proportion of the union (U) of both distributions. All values were converted to percentages by multiplying by 100.

$$ \text{TNO}(\%) = 100 \times \frac{T_{\text{recent}} \ominus T_{\text{baseline}}}{T_{\text{recent}} \cup T_{\text{baseline}}} = 100 \times \left(1 - \frac{T_{\text{recent}} \cap T_{\text{baseline}}}{T_{\text{recent}} \cup T_{\text{baseline}}}\right) \quad (1) $$

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<sub>o</sub>) 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.

$$ \text{TNO}_{o}(\%) = 100 \times \frac{T_{\text{random } a} \ominus T_{\text{random } b}}{T_{\text{random } a} \cup T_{\text{random } b}} \quad (2) $$

To calculate standardized thermal non-overlap (TNO<sub>o</sub>), we subtracted TNO<sub>o</sub> from TNO thereby setting the null expectation to zero.

$$ \text{TNO}_{s}(\%) = \text{TNO} - \text{TNO}_{o} \quad (3) $$

In this case, if the temperature distributions in the recent and baseline time periods were identical, the TNO<sub>o</sub> 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<sub>o</sub>. A comparison between raw values of TNO and TNO<sub>o</sub> 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<sup>68</sup> 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<sub>o</sub> (we considered a range of values from 5 to 10,000).
To assess the effect of seasonal habitat restrictions ($S_{\text{limit}}$) and volumetric habitat restrictions ($V_{\text{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³ so that habitat restrictions could be more readily compared across lakes. Comparing a $V_{\text{limit}}$ value of 0.8 across lakes of different sizes assumes that a habitat restriction of 2 m³ in a 10 m³ lake would be comparable to a 20 m³ habitat delineation in a 200 m³ lake. The actual size of the seasonal habitat restrictions for each lake in units of days were calculated using the value of $S_{\text{limit}}$ as follows:

$$S = (\text{doy}_{\text{max}} - \text{doy}_{\text{min}}) (1 - V_{\text{limit}})$$

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_{\text{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_{\text{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³ were calculated using the value of $V_{\text{limit}}$ as follows:

$$V = (\text{volume}) \times (1 - V_{\text{limit}})$$

where V is the volumetric habitat restriction in units of m³, volume is the lake's total volume and $V_{\text{limit}}$ is the volumetric habitat restriction value scaled from 0 to 0.95. For example, if a lake with a volume of 100 m³ had a $V_{\text{limit}}$ value of 0.75, we randomly selected temperature values from time periods a and b which were within four 25 m³ ($100 m^3 \times (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_{\text{limit}}$ and $V_{\text{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_{\text{limit}}$=0, $V_{\text{limit}}$=0.8), 'Species 2' ($S_{\text{limit}}$=0.8, $V_{\text{limit}}$=0) and 'Species 3' ($S_{\text{limit}}$=0.8, $V_{\text{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_{\text{limit}}$=0.74, $V_{\text{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⁴¹. Of the available lake characteristics, we used latitude and mean depth because they were most strongly correlated to $TNO$ 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⁶⁸. 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²⁴. 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_{\text{limit}}$), depth habitat restrictions ($V_{\text{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$, 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⁶⁹.

**Reporting Summary.** Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

The remote sensing lake temperature data used can be found at https://doi.org/10.5067/GHAAO-4BC02 and through the Earthdata website (https://earthdata.nasa.gov/). The in situ lake temperature data used are available through the Environmental Data Initiative (EDI) data portal (https://portalega repository.nis/home.jsp#) at https://doi.org/10.6073/pasta/f03d7d682ae2f467642e4260686ea15. Geospatial and morphometric data for each lake are available from the HydroLAKES database at https://doi.org/10.1038/ncomms13603 and can be found at http://www.hydrosheds.org. The country vector data shown in Fig. 3 were sourced from Natural Earth Data (https://www.naturalearthdata.com) under a public domain licence that allows full and unrestricted use. The lake-specific panels in Figs. 2 and 4 have been repeated for all 139 lakes and are available at https://doi.org/10.5281/zenodo.4688294.

## Code availability

All code is available under the identifier https://doi.org/10.5281/zenodo.4688292.

## References

57. Livingstone, D. M. Impact of secular climate change on the thermal structure of a large temperate central European lake. *Climatic Change* 57, 205–225 (2003).

58. Coats, R., Perez-Losada, J., Schladow, G., Richards, R. & Goldman, C. The warming of Lake Tahoe. *Climatic Change* 76, 121–148 (2006).

59. Gaiser, E. E., Deyrup, N. D., Bachmann, R. W., Battoe, L. E. & Swain, H. M. Multidecadal climate oscillations detected in a transparency record from a subtropical Florida lake. *Linnol. Oceanogr.* 54, 2228–2232 (2009).

60. Pilla, R.M. et al. Global data set of long-term summertime vertical temperature profiles in 153 lakes. *Environ. Data Initiat.* https://doi.org/10.6073/pasta/f03d7d682ae2f467642e4260686ea15 (2021).

61. Rimet, F. et al. The Observatory on LAKES (OLA) database: sixty years of environmental data accessible to the public. *J. Limnol.* 79, 2 (2020).

62. MacIntyre, S. *Water Temperatures from Toolkit Lake in Summer 2004* (Environmental Data Initiative, 2015); https://doi.org/10.6073/pasta/3fdd4c2d590d5485770e031d7e8678d8 (2006).

63. Kraemer, B. *Century-Long Temperature Record from Lake Tanganyika* (Environmental Data Initiative, 2015); https://doi.org/10.6073/pasta/59946ad9f1feffd541d60654490b104

64. N. Lead PI et al. *North Temperate Lakes LTER: Physical Limnology of Primary Study Lakes 1981 - Current* (Environmental Data Initiative, 2010); https://doi.org/10.6073/pasta/c120b223f80c63982457a2e1e76f6038
65. NCEI. GHRSST Level 4 AVHRR_OI Global Blended Sea Surface Temperature Analysis (GDS v.2) from NCEI. v.2.0 (PO.DAAC, 2016); https://doi.org/10.5067/GHAAO-4BC02

66. Batt, R. Free-water lake metabolism: addressing noisy time series with a Kalman filter. *Limnol. Oceanogr.-Meth.* **10**, 20–30 (2012).

67. Johansson, H., Brolin, A. A. & Håkanson, L. New approaches to the modelling of lake basin morphometry. *Environ. Model. Assess.* **12**, 213–228 (2007).

68. Elith, J., Leathwick, J. & Hastie, T. A working guide to boosted regression trees. *J. Anim. Ecol.* **77**, 802–813 (2008).

69. R Core Team. *R: A Language and Environment for Statistical Computing* (R Foundation for Statistical Computing, 2017).

## Acknowledgements

This work was initiated through the Global Lake Ecological Observatory Network (GLEON) and benefitted from continued participation and travel support from GLEON. We thank C. Munteanu, K. Kakouei, M. Thayne and U. Scharfenberger for their constructive discussions about the manuscript. Funding and data in support of this work also came from the following sources: Belarus Republican Foundation for Fundamental Research; Leibniz Institute for Freshwater Ecology and Inland Fisheries Long-Term Research (Germany); 2017–2018 Belmont Forum and BiodivERsA joint call for research proposals under the BiodivScen ERA-Net COFUND programme; and with the funding organizations German Science Foundation (grant no. AD 91/22-1); Universidad del Valle de Guatemala; Oklahoma Department of Wildlife Conservation (through the Sport Fish Restoration Program, Grant F-61-R; United States); Archbold Biological Station; Oklahoma Water Resources Board; Grand River Dam Authority; United States Army Corps of Engineers; City of Tulsa; Ministry of Business, Innovation and Employment (grant no. UOW X1503; New Zealand); Natural Environment Research Council of the United Kingdom (grant no. NE/R016429/1); Lacawac Sanctuary and Biological Field Station; Ontario Ministry of the Environment; Natural Sciences and Engineering Research Council of Canada; Canada Foundation for Innovation; Canada Research Chairs; Province of Saskatchewan, Canada; University of Regina, Canada; Queen's University Belfast, United Kingdom; California Air Resources Board; United States National Aeronautics and Space Administration; United States National Park Service; European Ministry of Higher Education and Research (grant no. 671 6.1387.2017); Environmental Agency of Verona; United States National Science Foundation (grant nos. DEB-1754276, DEB-1242626 and DEB-1950170); Gordon and Betty Moore Foundation; Mellon Foundation; University of Washington; Tyrolean Alps Long-Term Sociological Ecological and Research (LTSER, Austria); Waikato Regional Council; Bay of Plenty Regional Council; the Observatory of Lakes (OLA); SILA; CISALB; CIPEL (SOERE-OLA, [downloaded in January 2014]); The Norwegian Water Resources and Energy Directorate (NVE), by courtesy of Å. S. Kvambekk; Russian Science Foundation

(grant no. 20-64-46003); Russian Ministry of Higher Education and Research (grant nos. FZZE-2020-0026 and FZZE-2020-0023); Lake Baikal Foundation for Support of Applied Ecological Research and Development (https://baikalfoundation.ru/project/tochka-1/); British Antarctic Survey; International Long-Term Ecological Research Network; United States Environmental Protection Agency; Lake Tahoe Environmental Research Center; Yigal Allon Kinneret Limnological Laboratory; City of Zurich Water Supply; Office of Waste, Water, Energy and Air, Canton of Zürich; Shiga Prefectural Fishery Experiment Station; Experimental Lakes Area of Canada; University of Nevada, Reno; Belgian Federal Science Policy Office (BELSPO project EAGLES CD/AR/02A); National Capital Authority (ACT, Australia); Swedish Environmental Protection Agency; Swedish Infrastructure for Ecosystem Sciences financed by the Swedish Research Council (grant no. 2017-00635); International Commission for the Protection of Italian-Swiss Waters (CIPAIS); International Commission for the Protection of Lake Constance; Max Planck Institute for Limnology; and the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant no. 722518 (MANTEL).

## Author contributions

B.M.K. conceived of the work, led the study, completed the analysis and wrote the manuscript with input from all co-authors. All authors edited and revised the manuscript. R.A., P.R.L., R.M.P. and R.I.W. provided additional feedback, edits and revisions. Data contributions were made by O.A., S.B., W.C.-M., S.P.D., M.T.D., E.E.G., K.D.H., D.O.H., S.N.H., K.D.J., W.K., L.B.K., B.M.K., P.R.L., F.L., M.S.L., S.C.M., D.C.M.-N., D.C.P., R.M.P., D.C.R., M.R., A.M.P., J.A.R., S.S., N.S., E.A.S., D.S., M.A.T., W.T., P.V., R.S., G.A.W. and R.A. R.M.P. and B.M.K. compiled, organized and quality checked the data across all lakes.

## Competing interests

The authors declare no competing interests.

## Additional information

**Extended data** is available for this paper at https://doi.org/10.1038/s41558-021-01060-3.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41558-021-01060-3.

**Correspondence and requests for materials** should be addressed to B.M.K.

**Peer review information** *Nature Climate Change* thanks Gretchen Hansen, Conor Waldock and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.
<figure><img src="image_27.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_30.png" /></figure>

Extended Data Fig. 1 | See next page for caption.
**Extended Data 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.
<figure><img src="image_31.png" /></figure>

**Extended Data 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-I) with the most important predictors in the upper left (a) and the least important in the bottom right (I). Note that the scales on the x- and y-axes vary from plot to plot.
<figure><img src="image_32.png" /></figure>

**Extended Data Fig. 3** | Relative importance of each of the 12 predictors in the boosted regression trees. The relative importance of each predictor variable in the model is a function of the frequency with which it was included in the BRT's individual regression trees and the improvement to the model that resulted from its inclusion. Relative importance values are scaled to sum to 100%. The individual colored points represent 100 different BRTs fit to random subsets (70%) of the 139 lakes included in this analysis and they are jittered to aid their visualization. The colored violins are the kernel density estimations of the distribution and the black boxes show the first, second and third quartiles of the distributions.
<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /><figcaption>Extended Data Fig. 4 | See next page for caption.</figcaption></figure>
**Extended Data Fig. 4 | Comparison of thermal non-overlap and whole-lake warming rates across all 139 lakes.** The density plot (a) shows the distribution of volume-weighted whole-lake warming rates across the 139 lakes included in the study. Thermal non-overlap and the absolute value of whole-lake warming rates were not strongly correlated (r = 0.10, p =0.081, Kendall's rank correlation; b).
<figure><img src="image_35.png" /></figure>

**Extended Data Fig. 5 | Comparison of raw thermal non-overlap values with null non-overlap values prior to removing the effects of time series characteristics.** Raw thermal non-overlap values were calculated by comparing the first and second parts of each lake's time series whereas the null non-overlap values were calculated by comparing two groups of randomly selected years. Values shown here are for the case with no seasonal or depth habitat restrictions. The non-overlap values in the main text are 'standardized' by subtracting the null non-overlap values from the raw non-overlap values shown here.
<table><thead><tr><th>Summary statistic</th><th>Value</th></tr></thead><tbody><tr><td>Training data correlation between predicted and observed values where correlations were averaged across all BRTs</td><td>Kendall = 0.96<br>Spearman = 0.99</td></tr><tr><td>Training data correlation between predicted and observed values where predicted values were averaged across all BRTs</td><td>Kendall = 0.98<br>Spearman = 0.99</td></tr><tr><td>Test data correlation between predicted and observed values where correlations were averaged across all BRTs</td><td>Kendall = 0.56<br>Spearman = 0.75</td></tr><tr><td>Test data correlation between predicted and observed values where predicted values were averaged across all BRTs</td><td>Kendall = 0.61<br>Spearman = 0.81</td></tr><tr><td>Overall correlation between predicted and observed values where correlations were averaged across all BRTs and individual correlations were calculated from pooled training and test datasets</td><td>Kendall = 0.92<br>Spearman = 0.93</td></tr><tr><td>Overall correlation between predicted and observed values where predicted values were averaged across all BRTs from pooled training and test datasets</td><td>Kendall = 0.90<br>Spearman = 0.99</td></tr><tr><td>Standard deviation (SD) of standardized thermal non-overlap values (TNOs) prior to cross validation</td><td>16.24 % non-overlap</td></tr><tr><td>Training data SD of model residuals where SDs were averaged across all BRTs</td><td>0.99 % non-overlap</td></tr><tr><td>Training data SD of model residuals where predicted values were averaged across all BRTs</td><td>0.59 % non-overlap</td></tr><tr><td>Test data SD of model residuals (PRESS) where SDs were averaged across all BRTs</td><td>10.44 % non-overlap</td></tr><tr><td>Test data SD of model residuals (PRESS) where predicted values were averaged across all BRTs</td><td>9.36 % non-overlap</td></tr><tr><td>Overall SD between predicted and observed values where SDs were averaged across all BRTs and individual SDs were calculated from pooled training and test datasets</td><td>5.84 % non-overlap</td></tr><tr><td>Overall SD between predicted and observed values where predicted values were averaged across all BRTs from pooled training and test datasets</td><td>2.93 % non-overlap</td></tr><tr><td>Median absolute deviation (MAD) of standardized thermal non-overlap values (TNOs) prior to cross validation</td><td>11.34 % non-overlap</td></tr></tbody></table>

**Extended Data Fig. 6** | Model performance statistics from the cross-validation for the boosted regression trees analysis. Summary statistics are described in the 'Summary statistic' column and the values of each statistic are reported in the adjacent 'value column'.
<figure><img src="image_36.png" /></figure>

**Extended Data Fig. 7** | Analysis of model residuals calculated as the difference between observed values and the average predicted values across all 100 BRTs. The distribution of thermal habitat change (% non-overlap) values compared to the distribution of model residuals (a). Model residuals show minimal patterning when compared to each predictor in the BRT analysis (b:m). Colors represent the counts of observations in 2-dimensional bins where the color scale has been log₁₀-transformed.
# Reporting Summary

Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see our Editorial Policies and the Editorial Policy Checklist.

## Statistics

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

<figure><img src="image_37.png" /></figure>

Our web collection on [statistics for biologists](#) contains articles on many of the points above.

## 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.
# 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.

<figure><img src="image_38.png" /></figure>

Life sciences

<figure><img src="image_39.png" /></figure>

Behavioural & social sciences

<figure><img src="image_40.png" /></figure>

Ecological, evolutionary & environmental sciences

For a reference copy of the document with all sections, see [nature.com/documents/nr-reporting-summary-flat.pdf](https://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.

**Blinding**
Blinding was not relevant for this observational study based on long-term monitoring data which does not include experiments.

Did the study involve field work?

<figure><img src="image_41.png" /></figure>

Yes

<figure><img src="image_42.png" /></figure>

No

# Field work, collection and transport

**Field conditions**
Lake water temperature measurements were collected across 139 lakes under a wide variety of field conditions over several decades.

**Location**
All sampling locations can be found in the supplementary material for the manuscript.

**Access & import/export**
All temperature data were recorded digitally or in writing upon measurement and no physical samples were imported or exported for the purpose of this study.

**Disturbance**
The measurement of lake water temperature is generally not an intrusive process.

# Reporting for specific materials, systems and methods

We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.
Materials & experimental systems

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☑</td><td>Antibodies</td></tr><tr><td>☑</td><td>Eukaryotic cell lines</td></tr><tr><td>☑</td><td>Palaeontology and archaeology</td></tr><tr><td>☑</td><td>Animals and other organisms</td></tr><tr><td>☑</td><td>Human research participants</td></tr><tr><td>☑</td><td>Clinical data</td></tr><tr><td>☑</td><td>Dual use research of concern</td></tr></tbody></table>

Methods

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☑</td><td>ChIP-seq</td></tr><tr><td>☑</td><td>Flow cytometry</td></tr><tr><td>☑</td><td>MRI-based neuroimaging</td></tr></tbody></table>