# Widespread deoxygenation of temperate lakes

Article

Accepted Version

Jane, S. F., Hansen, G. J. A., Kraemer, B. M., Leavitt, P. R., Mincer, J. L., North, R. L., Pilla, R. M., Stetler, J. T., Williamson, C. E., Woolway, R. I. ORCID: https://orcid.org/0000-0003-0498-7968, Arvola, L., Chandra, S., DeGasperi, C. L., Diemer, L., Dunalska, J., Erina, O., Flaim, G., Grossart, H.-P., Hambright, K. D., Hein, C., Hejzlar, J., Janus, L. L., Jenny, J.-P., Jones, J. R., Knoll, L. B., Leoni, B., Mackay, E., Matsuzaki, S.-i. S., McBride, C., Müller-Navarra, D. C., Paterson, A. M., Pierson, D., Rogora, M., Rusak, J. A., Sadro, S., Saulnier-Talbot, E., Schmid, M., Sommaruga, R., Thiery, W., Verburg, P., Weathers, K. C., Weyhenmeyer, G. A., Yokota, K. and Rose, K. C. (2021) Widespread deoxygenation of temperate lakes. Nature, 594. pp. 66-70. ISSN 0028-0836 doi: https://doi.org/10.1038/s41586-021-03550-y Available at https://centaur.reading.ac.uk/100084/

It is advisable to refer to the publisher's version if you intend to cite from the work. See Guidance on citing.

To link to this article DOI: http://dx.doi.org/10.1038/s41586-021-03550-y

Publisher: Nature Publishing Group
All outputs in CentAUR are protected by Intellectual Property Rights law,
including copyright law. Copyright and IPR is retained by the creators or other
copyright holders. Terms and conditions for use of this material are defined in
the [End User Agreement](#).

www.reading.ac.uk/centaur

CentAUR

Central Archive at the University of Reading

Reading's research outputs online
1 **Title:**

2 Widespread deoxygenation of temperate lakes

3

4 **Authors:**

5 Stephen F. Jane¹,², Gretchen J. A. Hansen³, Benjamin M. Kraemer⁴, Peter R. Leavitt⁵,⁶, Joshua L.

6 Mincer¹, Rebecca L. North⁷, Rachel M. Pilla⁸, Jonathan T. Stetler¹, Craig E. Williamson⁸, R.

7 Iestyn Woolway⁹, Lauri Arvola¹⁰, Sudeep Chandra¹¹, Curtis L. DeGasperi¹², Laura Diemer¹³,

8 Julita Dunalska¹⁴, Oxana Erina¹⁵, Giovanna Flaim¹⁶, Hans-Peter Grossart¹⁷, ¹⁸, K. David

9 Hambright¹⁹, Catherine Hein²⁰, Josef Hejzlar²¹, Lorraine L. Janus²², Jean-Philippe Jenny²³, John

10 R. Jones⁷, Lesley B. Knoll²⁴, Barbara Leoni²⁵, Eleanor Mackay²⁶, Shin-Ichiro S. Matsuzaki²⁷,

11 Chris McBride²⁸, Dörthe C. Müller-Navarra²⁹, Andrew M. Paterson³⁰, Don Pierson², Michela

12 Rogora³¹, James A. Rusak²⁹, Steven Sadro³², Emilie Saulnier-Talbot³³, Martin Schmid³⁴, Ruben

13 Sommaruga³⁵, Wim Thiery³⁶, ³⁷, Piet Verburg³⁸, Kathleen C. Weathers³⁹, Gesa A.

14 Weyhenmeyer², Kiyoko Yokota⁴⁰, and Kevin C. Rose¹

15

16 1. Department of Biological Sciences, Rensselaer Polytechnic Institute, Troy, NY, USA

17 2. Department of Ecology and Genetics/Limnology, Uppsala University, Uppsala, Sweden

18 3. Department of Fisheries, Wildlife and Conservation Biology, University of Minnesota,
19 St. Paul, MN, USA

20 4. Department of Ecosystem Research, IGB Leibniz institute for freshwater ecology and
21 inland fisheries, Berlin, Germany
<table>
  <tr>
    <td>22</td>
    <td>5. Institute of Environmental Change and Society, University of Regina, Regina, Saskatchewan, Canada</td>
  </tr>
  <tr>
    <td>23</td>
    <td>6. Institute for Global Food Security, Queen's University Belfast, Belfast, County Antrim, United Kingdom</td>
  </tr>
  <tr>
    <td>24</td>
    <td>7. School of Natural Resources, University of Missouri, Columbia, MO, USA</td>
  </tr>
  <tr>
    <td>25</td>
    <td>8. Department of Biology, Miami University, Oxford, OH, USA</td>
  </tr>
  <tr>
    <td>26</td>
    <td>9. Centre for Freshwater and Environmental Studies, Dundalk Institute of Technology, Dundalk, Ireland</td>
  </tr>
  <tr>
    <td>27</td>
    <td>10. Lammi Biological Station, University of Helsinki, Pääjärventie 320, Lammi, Finland</td>
  </tr>
  <tr>
    <td>28</td>
    <td>11. Biology Department and Global Water Center, University of Nevada, Reno, NV, USA</td>
  </tr>
  <tr>
    <td>29</td>
    <td>12. King County Water and Land Resources Division, Seattle, WA, USA</td>
  </tr>
  <tr>
    <td>30</td>
    <td>13. FB Environmental Associates, Portsmouth, NH, USA</td>
  </tr>
  <tr>
    <td>31</td>
    <td>14. Department of Water Protection Engineering, University of Warmia and Mazury in Olsztyn, Prawocheńskiego str. 1, Olsztyn, Poland</td>
  </tr>
  <tr>
    <td>32</td>
    <td>15. Department of Hydrology, Lomonosov Moscow State University, Leninskiye Gory 1, Moscow, Russia</td>
  </tr>
  <tr>
    <td>33</td>
    <td>16. Department of Sustainable Agro-ecosystems and Bioresources, Research and Innovation Centre, Fondazione Edmund Mach, San Michele all'Adige, Italy</td>
  </tr>
  <tr>
    <td>34</td>
    <td>17. Department of Experimental Limnology, Leibniz-Institute of Freshwater Ecology and Inland Fisheries, Alte Fischerhuette 2, Stechlin, Germany</td>
  </tr>
  <tr>
    <td>35</td>
    <td>18. Department of Water Protection Engineering, University of Warmia and Mazury in Olsztyn, Prawocheńskiego str. 1, Olsztyn, Poland</td>
  </tr>
  <tr>
    <td>36</td>
    <td>19. Department of Hydrology, Lomonosov Moscow State University, Leninskiye Gory 1, Moscow, Russia</td>
  </tr>
  <tr>
    <td>37</td>
    <td>20. Department of Water Protection Engineering, University of Warmia and Mazury in Olsztyn, Prawocheńskiego str. 1, Olsztyn, Poland</td>
  </tr>
  <tr>
    <td>38</td>
    <td>21. Department of Hydrology, Lomonosov Moscow State University, Leninskiye Gory 1, Moscow, Russia</td>
  </tr>
  <tr>
    <td>39</td>
    <td>22. Department of Water Protection Engineering, University of Warmia and Mazury in Olsztyn, Prawocheńskiego str. 1, Olsztyn, Poland</td>
  </tr>
  <tr>
    <td>40</td>
    <td>23. Department of Experimental Limnology, Leibniz-Institute of Freshwater Ecology and Inland Fisheries, Alte Fischerhuette 2, Stechlin, Germany</td>
  </tr>
  <tr>
    <td>41</td>
    <td>24. Department of Water Protection Engineering, University of Warmia and Mazury in Olsztyn, Prawocheńskiego str. 1, Olsztyn, Poland</td>
  </tr>
</table>
42 18. Institute of Biochemistry and Biology, Potsdam University, Maulbeerallee 2, Potsdam,
43 Germany

44 19. Plankton Ecology and Limnology Laboratory, Geographical Ecology Group, and
45 Program in Ecology and Evolutionary Biology, Department of Biology, The University
46 of Oklahoma, Norman, Oklahoma, USA

47 20. Wisconsin Department of Natural Resources, Madison, WI, USA

48 21. Institute of Hydrobiology, Biology Centre CAS, Na Sádkách 7, České Budějovice, Czech
49 Republic

50 22. Bureau of Water Supply, New York City Department of Environmental Protection, 465
51 Columbus Ave., Valhalla, NY, USA

52 23. CARRTEL Limnology Center, Institut National de la Recherche Agronomique (INRA),
53 Université Savoie Mont Blanc, 73000 Chambéry, France

54 24. Itasca Biological Station and Laboratories, University of Minnesota, Lake Itasca, MN,
55 USA

56 25. Department of Earth and Environmental Sciences, University of Milan-Bicocca, Piazza
della Scienza 1, Milan, Italy

58 26. Lake Ecosystems Group, Centre for Ecology & Hydrology, Bailrigg, Lancaster, United
59 Kingdom

60 27. Center for Environmental Biology and Ecosystem Studies, National Institute for
61 Environmental Studies, 16-2 Tsukuba, Ibaraki, Japan

62 28. Environmental Research Institute, Hamilton, New Zealand
63 29. Department of Biology, University of Hamburg, Hamburg, Germany

64 30. Ontario Ministry of the Environment, Conservation and Parks, Dorset Environmental
65 Science Centre, Dorset, ON, Canada

66 31. CNR Water Research Institute, L. go Tonolli 50, Verbania Pallanza, Italy

67 32. Department of Environmental Science and Policy, University of California, One Shields
68 Ave., Davis, CA, USA

69 33. Centre d'études nordiques (CEN) and Faculty of Sciences and Engineering, Université
70 Laval, Québec, Canada

71 34. Eawag, Swiss Federal Institute of Aquatic Science and Technology, Surface Waters –
72 Research and Management, Kastanienbaum, Switzerland

73 35. Department of Ecology, University of Innsbruck, Innsbruck, Austria

74 36. Vrije Universiteit Brussel, Department of Hydrology and Hydraulic Engineering,
75 Pleinlaan 2, Brussels, Belgium

76 37. ETH Zurich, Institute for Atmospheric and Climate Science, Universitaetstrasse 16,
77 Zurich, Switzerland

78 38. National Institute of Water and Atmospheric Research Ltd (NIWA), Gate 10 Silverdale
79 Rd., Hillcrest, Hamilton, New Zealand

80 39. Cary Institute of Ecosystem Studies, Box AB, Millbrook, New York, USA

81 40. Biology Department, State University of New York College at Oneonta (SUNY
82 Oneonta), Oneonta, New York, USA
83 **Summary paragraph:**

The concentration of dissolved oxygen in aquatic systems helps regulate biodiveristy¹, ², nutrient biogeochemistry³, greenhouse gas emissions⁴, and drinking water quality⁵. The long-term declines in dissolved oxygen concentrations in coastal and ocean waters have been linked to climate warming and human activity⁶, ⁷, but little is known about changes in dissolved oxygen concentrations in lakes. While dissolved oxygen solubility decreases with increasing water temperatures, long-term lake trajectories are not necessarily predictable. Oxygen losses in warming lakes may be amplified by enhanced decomposition and stronger thermal stratification⁸, ⁹ or they may increase as a result of enhanced primary production¹⁰. Here we analyse 45,148 dissolved oxygen and temperature profiles from 393 temperate lakes spanning 1941-2017. We find that a decline in dissolved oxygen is widespread in surface and deep-water habitats. The decline in surface waters is primarily associated with reduced solubility under warmer water temperatures, although surface dissolved oxygen increased in a subset of highly-productive warming lakes, likely due to increasing phytoplankton production. In contrast, the decline in deep waters is associated with stronger thermal stratification and water clarity losses, but not with changes in gas solubility. Our results suggest that climate change and declining water clarity have altered the physical and chemical environment of lakes. Freshwater dissolved oxygen losses are 2.5-10 times greater than observed in the world's oceans⁶, ⁷ and could threaten essential lake ecosystem services², ³, ⁵, ¹¹.

102

103

104
105

106

107

108

109

110

111

112

113

The concentration of dissolved oxygen (DO) in aquatic systems influences biodiversity¹,
², nutrient biogeochemistry³, greenhouse gas emissions⁴, drinking water quality⁵, and, ultimately,
human health¹². Many aquatic species require well-oxygenated habitat¹¹, ¹³ and cool water to
survive warm summers², ¹¹. Loss of deep-water DO degrades water quality by promoting the
release of accumulated nutrients from sediments into water¹, ³, which can increase phytoplankton
biomass. This process can also facilitate harmful algal blooms⁵, which can compromise water
supplies and harm human health¹². Despite clear evidence of large-scale deoxygenation in ocean
waters⁶, ⁷, there are no systematic large-scale studies of this phenomenon in lakes³.

114

115

116

117

118

119

120

121

122

DO concentrations should decline with increasing water temperature due to reduced gas solubility. However, other mechanisms can alter DO, potentially amplifying or counteracting losses predicted from solubility changes alone. For example, rates of heterotrophic respiration increase with temperature faster than primary production⁹, and surface-temperature warming can increase the strength and duration of thermal stratification, reducing water circulation, and preventing deep-water DO replenishment⁸, ¹⁴, ¹⁵. Studies of individual lakes demonstrate deep-water DO concentrations can decrease with lake warming³, ⁸, ¹⁵, ¹⁶, reducing access to cold-water habitat essential to many organisms¹¹. However, given the many feedbacks and processes regulating DO, overall trajectories currently defy a priori prediction.

123

124

125

126

127

We addressed this critical issue by compiling and analyzing an extensive database of lake
temperature and DO profiles to characterize widespread and long-term changes in DO
concentration and its causes. We used data from 393 temperate lake and reservoir basins, each
with a minimum of 15 years of observation (median: 24 years), and report population medians
from long-term surface- (epilimnion) and deep-water (hypolimnion) trends in temperature, DO
concentration, and DO saturation during the late summer period when seasonal DO depletion is expected to be pronounced¹⁷. Our analyses revealed that lake DO concentrations have declined in both surface and deep waters from 1980 to 2017 by 0.45 and 0.42 mg L⁻¹, respectively (Fig. 1). These rates represent losses of 5.1 and 20.2% for surface and deep waters, respectively, and were substantially greater than those observed for the oceans, where total water-column DO has declined about 2% since 1960⁶.

While deep-water temperatures have been virtually stable since observations began (Fig. 1a; -0.01°C decade⁻¹), both deep-water DO concentration and percent saturation declined through time (-0.12 mg L⁻¹ decade⁻¹ and -1.2% decade⁻¹; respectively, Fig. 1b, c). Declines were unrelated to solubility as predicted changes based on solubility (slight increase of 0.01 mg L⁻¹) were negligible compared with observed losses (median -0.23 mg L⁻¹ based on last five years relative to first five years of each time series, Fig. 2b). Declining DO, despite essentially unchanging solubility, implies deep-water habitats have become increasingly inhospitable for organisms with aerobic metabolism, including fishes. We quantified potential impacts of such declines on habitat availability by calculating trends in TDO₃, the minimum water column temperature where DO was at least 3 mg L⁻¹. This metric was developed to quantify oxy-thermal habitats for cold-water fisheries¹¹. In lakes where DO was below 3 mg L⁻¹ anywhere in the water column at least once in the time series (n = 369), TDO₃ increased by 0.17°C decade⁻¹, with 68.0% of lakes having positive trends and declining habitat for many cold-water species.

In contrast to trends observed for deep waters, variation in surface-water DO concentrations was well explained by changes in gas solubility. Consistent with other global-scale lake studies¹⁸, median air temperatures warmed at 0.30°C decade⁻¹ and median lake surface waters warmed at 0.39°C decade⁻¹. Additionally, median wind speed and precipitation declined
(trends of -0.04 m s⁻¹ decade⁻¹ and -4.23 mm decade⁻¹, respectively), while shortwave radiation increased (1.88 W m² decade⁻¹; Table S1). Surface-water temperature increases were best explained by spring and summer air temperature increases and by summer wind speed declines (Table S2). Surface-water DO concentrations declined at -0.11 mg L⁻¹ decade⁻¹ (Fig. 1b).

Comparing the last five years relative to first five years of each time series revealed that the median change predicted due to solubility loss was ~63% of the median observed decline in DO concentration, with solubility-predicted loss of 0.12 versus observed losses of 0.19 mg L⁻¹ (Fig. 2a).

Despite a strong influence of water temperature on DO concentration in surface-waters, there was substantial variability among lakes (Fig. 2a), and a large subset of lakes exhibited increases in both water temperature and DO concentration (n=87; Fig. 3d). Analysis of the interaction between DO concentration, surface temperature, and water clarity (measured as Secchi depth, a proxy for trophic status¹⁹) showed that DO concentration generally decreased with increasing temperature. However, in lakes with low water clarity (< 2 m), DO concentration increased when average mean summer surface-water temperatures exceeded ~24°C (Fig. 3c).

Similarly, in a subset of lakes with chlorophyll data (a proxy for phytoplankton biomass; n = 162), positive DO trends were observed when chlorophyll was high and surface temperatures exceeded ~25°C, (Fig. 3b; *P* < 0.001). Thus, we suggest that eutrophication and warming interact to increase surface-water DO concentration despite reduced gas solubility.

Lakes with increasing DO concentration in warming surface waters had significantly higher surface-water temperatures (Fig. 3a; *P* = 0.016) and their watersheds contained a significantly higher proportion of agriculture (*P* = 0.046) and developed land cover (*P* < 0.001) compared with other lakes. When developed land exceeded ~50% of a watershed and surface
water temperature exceeded ~25°C, the probability of a warming lake having an increasing DO trend was >50%. Combined, these analyses highlight a potential threshold above which water temperatures and lake productivity interact to elevate DO concentration in surface waters despite declining gas solubility. While we lack data on phytoplankton taxonomic composition, evidence indicates that phytoplankton blooms are increasing globally²⁰, in particular due to cyanobacteria²¹. High temperatures and elevated nutrient loading can promote surface cyanobacteria blooms whose photosynthesis leads to DO supersaturation, particularly in eutrophic lakes as temperatures exceed ~23-25°C¹⁰, ²¹. Consistent with this inferred mechanism, we note these same lakes exhibited consistently low deep-water DO concentration (median: 0.64 mg L⁻¹) relative to other lakes (median: 3.42 mg L⁻¹), as is expected when a large phytoplankton biomass sinks and is decomposed in deep-water habitats²². Deep water DO changes are described in more detail below.

Decadal-scale trends in DO were associated with non-linear changes in surface-water temperature (Fig. 2c-f; Fig. S1). For example, although surface-water temperatures generally increased from 1980 onwards, there was a period of accelerated increase during 1990-2000, with slower warming thereafter (Fig. 2c), consistent with the “warming hiatus” observed during 1998-2012²³. This trend occurs across the population of all lakes, as well as the subset of lakes sampled continuously throughout this period. Similarly, surface-water DO exhibited periodic deviations from an overarching trend of declining DO concentration (Fig. 2d), mainly due to the productive lakes exhibiting increasing DO levels in surface waters (Fig. 2d, blue line). Excluding these lakes, analysis of the remaining sites showed a consistent long-term decline in surface-water DO (Fig. 2d, red line). Deep-water temperatures exhibited a pronounced multi-decadal
oscillation since 1980 (Fig. 2e) as has been observed in some lakes previously²⁴, whereas deep-water DO concentration declined consistently through time (Fig. 2f).

While surface-water DO concentration changes were generally well predicted by solubility changes, deep-water DO changes were more strongly associated with changes in water clarity and water-column density differences (Figs. 4 and S2). For example, water clarity losses exceeding 1 m were associated with substantial reductions in deep-water DO saturation (Fig. S2). Mechanistically, increases in phytoplankton biomass or dissolved organic matter (DOM) reduce water clarity while increasing oxygen-consuming respiration¹⁹, ²², ²⁵. Increases in phytoplankton biomass and DOM are often caused by land use change and recovery from acid deposition, respectively²⁶. However, there was no overarching decline in water clarity across study lakes. Indeed, 51% of lakes had clarity increases and 49% had decreases, and only 39% of lakes exhibited both water clarity loss and DO saturation loss (Fig. 4a).

Deep-water DO decreased substantially in lakes where the water column density difference between surface and deep waters increased by more than ~0.5 kg m⁻³ (Fig 4b; Fig. S2b). Strong increases in the density difference indicate intensified stratification that reduces vertical mixing and replenishment of deep-water DO from the atmosphere, and may reduce nutrient upwelling to surface waters³, ¹⁵. Water column density differences increase due to water clarity losses as well as other factors that increase heat gain in near-surface waters, including climate warming²⁶ and atmospheric stilling²⁷. Increased water column density differences may also be associated with earlier onset of seasonal stratification and thus more time for oxygen consumption before the summer sampling period²². We found that changes in water-column density differences were best explained by changes in deep water temperature and climate characteristics (Fig. S3). Despite no overarching among-lake trend in water clarity or deep-water
temperature, stratification strength increased in 84% of lakes that stratified, with 61% of basins exhibiting both increased density difference and DO saturation loss (Fig 4b). Warming surface-water temperatures combined with unchanging deep-water temperatures (Fig. 1a) increases the density difference in lake water columns (median rate: 0.10 kg m⁻³ decade⁻¹). We observed unchanging deep-water DO in lakes where both clarity and stratification were unchanged (Fig. 4c, d). Therefore, we anticipate further DO losses in deep waters of lakes where water clarity continues to decline or thermal stratification intensifies, whether due to atmospheric warming, stilling, or both²⁶, ²⁷.

Despite a wide range of lake and catchment characteristics, the overall trend of temperate lake deoxygenation is clear, with climate changes and water clarity losses contributing to declines in lake DO concentration at rates ~2.5-10 times greater than those observed in the global oceans⁶, ⁷. We find deep-water lake habitats are especially threatened, and deep-water DO trends may portend future losses of cold-water and oxygen-sensitive species², increased internal nutrient loading which exacerbates eutrophication³ and the formation of harmful algal blooms⁵, and potentially increased outgassing of stored methane⁴. While already rapid, future losses in lake DO may accelerate due to continued anthropogenic modifications of the environment, including eutrophication²², salinization²⁸, and hydrological management²⁸. While many lakes have undergone active management to reduce nutrient loads, in part to mitigate phytoplankton growth and deep-water oxygen loss²⁸, our findings suggest such actions will likely require more rigorous efforts in the future to counter the effects of climate and land use change.
240 **References:**

241 1. Wetzel, R. G. 2001. Chap. 9 oxygen. In: *Limnology*, 3rd edn (ed Wetzel R. G.), pp 151-168, Academic Press, San Diego.

243 2. Schindler, D. Warmer climate squeezes aquatic predators out of their preferred habitat. *Proc. Natl. Acad. Sci.*, **114**, 9764-9765 (2017).

245 3. North, R. P., North, R. L., Livingstone, D. M., Köster, O., & Kipfer, R. Long-term changes in hypoxia and soluble reactive phosphorus in the hypolimnion of a large temperate lake: consequences of a climate regime shift. *Glob. Change Biol.*, **20**, 811-823 (2014).

248 4. Fernández, J. E., Peeters, F., & Hofmann, H. Importance of the autumn overturn and anoxic conditions in the hypolimnion for the annual methane emissions from a temperate lake. *Environ. Sci. Technol.*, **48**, 7297-7304 (2014).

251 5. Michalak, A. M., et al. Record-setting algal bloom in Lake Erie caused by agricultural and meteorological trends consistent with expected future conditions. *Proc. Natl. Acad. Sci.*, **110**, 6448-6452 (2013).

254 6. Schmidtko, S., Stramma, L., & Visbeck, M. Decline in global oceanic oxygen content during the past five decades. *Nature*, **542**, 335-339 (2017).

256 7. Breitburg, D., et al. Declining oxygen in the global ocean and coastal waters. *Science*, **359**, DOI: 10.1126/science.aam7240 (2018).

258 8. Jankowski, J., Livingstone, D. M., Bührer, H., Forster, R., & Niederhauser, P. Consequences of the 2003 European heat wave for lake temperature profiles, thermal stability, and
260

hypolimnetic oxygen depletion: Implications for a warmer world. *Limnol. Oceanogr.*, **51**,

261

815-819 (2006).

262

9. Yvon-Durocher, G., Jones, J. I., Trimmer, M., Woodward, G., & Montoya, J. M. Warming

263

alters the metabolic balance of ecosystems. *Philos. T. R. Soc. B*, **365**, 2117-2126 (2010).

264

10. Seki, H., Takahashi, Y., Hara, Y., & Ichimura, S. Dynamics of dissolved oxygen during algal

265

bloom in Lake Kasumigaura, Japan. *Water Res.*, **14**, 179-183 (1980).

266

11. Jacobson, P. C., Stefan, H. G., & Pereira, D. L. Coldwater fish oxythermal habitat in

267

Minnesota lakes: influence of total phosphorus, July air temperature, and relative depth. *Can.*

268

*J. Fish. Aquat. Sci.*, **67**, 2002-2013 (2010).

269

12. Harke, M. J., et al. A review of the global ecology, genomics, and biogeography of the toxic

270

cyanobacterium, *Microcystis* spp. *Harmful Algae*, **54**, 4-20 (2016).

271

13. Vaquer-Sunyer, R., & Duarte, C. M. Thresholds of hypoxia for marine biodiversity. *Proc.*

272

*Natl. Acad. Sci.*, **105**, 15452-15457 (2008).

273

14. Woolway, R. I., & Merchant, C. J. Worldwide alteration of lake mixing regimes in response

274

to climate change. *Nat. Geosci.*, **12**, 271-276 (2019).

275

15. Livingstone, D. M. Impact of secular climate change on the thermal structure of a large

276

temperate central European lake. *Clim. Change*, **57**, 205-225 (2003).

277

16. Zhang, Y., et al. Dissolved oxygen stratification and response to thermal structure and long-

278

term climate change in a large and deep subtropical reservoir (Lake Qiandaohu, China).

279

*Water Res.*, **75**, 249-258 (2015).
17. Bouffard, D., Ackerman, J. D., & Boegman, L. Factors affecting the development and dynamics of hypoxia in a large shallow stratified lake: hourly to seasonal patterns. *Water Resour. Res.*, **49**, 2380-2394 (2013).

18. 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).

19. Nürnberg, G. K. Trophic state of clear and colored, soft- and hardwater lakes with special consideration of nutrients, anoxia, phytoplankton and fish. *Lake Reserv. Manag.*, **12**, 432-447 (1996).

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

21. Kosten, S., et al. Warmer climates boost cyanobacterial dominance in shallow lakes. *Glob. Change Biol.*, **18**, 118-126 (2012).

22. Müller, B., Bryant, L. D., Matzinger, A., & Wüest, A. Hypolimnetic oxygen depletion in eutrophic lakes. *Environ. Sci. Technol.*, **46**, 9964-9971 (2012).

23. Winslow, L. A., Leach, T. A., & Rose, K. C. Global lake response to the recent warming hiatus. *Environ. Res. Lett.*, **13**, 054005 (2018).

24. Livingstone, D. M. An example of the simultaneous occurrence of climate-driven “sawtooth” deep-water warming/cooling episodes in several Swiss lakes, *Verh. Int. Ver. Limnol.*, **26**, 822-828 (1997).

25. Williamson, C. E., et al. Ecological consequences of long-term browning in lakes. *Sci. Rep.*, **5**, DOI:10.1038/srep18666, (2015).
26. Rose, K. C., Winslow, L. A., Read, J. S., & Hansen, G. J. A. Climate-induced warming of lakes can be either amplified or suppressed by trends in water clarity. *Limnol. Oceanogr. Letters*, **1**, 44-53 (2016).

27. Woolway, R. I., Merchant, C. J., Van Den Hoek, J., Azorin-Molina, C., Nöges, P., Laas, A., Mackay, E. B., and Jones, I. D. Northern hemisphere atmospheric stilling accelerates lake thermal responses to a warming world. *Geophys. Res. Lett.*, **46**, 11983-11992 (2019).

28. Carpenter, S. R., Stanley, E. H., & Vander Zander, M. J. State of the world's freshwater ecosystems: physical, chemical, and biological changes. *Annu. Rev. Env. Resour.*, **36**, 75-99 (2011).

310

311
**Acknowledgments:** This manuscript benefited from conversations at meetings of the Global Lake Ecological Observatory Network (GLEON; supported by funding from US NSF grants 1137327 and 1702991). SFJ and KCR acknowledge support from US NSF grants 1638704, 1754265, and 1761805 and SFJ was supported by a US Fulbright Student grant to Uppsala University, Sweden. GH acknowledges the many employees of the Minnesota Department of Natural Resources, the Minnesota Pollution Control agency, and citizen volunteers for data collection and collation. PRL acknowledges support from a NSERC Discovery Grant, the Canada Research Chair Program, Canada Foundation for Innovation, the Province of Saskatchewan, University of Regina, and Queen's University Belfast. RLN & JJ acknowledge support from the Missouri Department of Natural Resources and the Missouri Agricultural Experiment Station and many students that collected and processed reservoir samples under the leadership of Daniel V. Obrecht and Anthony P. Thorpe. RMP and CEW acknowledge support from NSF grants 1754276 and 1950170, Miami University Eminent Scholar Fund, and the Lacawac Sanctuary and Biological Field Station for access to Lake Lacawac and use of research facilities. RIW acknowledges support from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 791812. SC acknowledges support of the Castle Lake Research Program through the University of Nevada and UC Davis via Charles R. Goldman. CLD acknowledges the King County Environmental Laboratory for the long-term monitoring data for Lake Washington and Lake Sammamish. JD acknowledges support from the University of Warmia and Mazury in Olsztyn (Grant under the Senate Committee for International Cooperation financing) and staff at Department of Water Protection Engineering for long-term data collection and analysis. OE acknowledges support from the Russian Scientific Foundation (grant 19-77-30004) for Mozhaysk Reservoir. GF
acknowledges support for long-term sampling of Lake Caldonazzo by the Fondazione Edmund Mach. HPG acknowledges funding for long-term sampling of Lake Stechlin by the Leibniz association and assistance by members of the IGB team. KDH acknowledges the Oklahoma Department of Wildlife Conservation, the Oklahoma Water Resources Board, the Grand River Dam Authority, the US Army Corps of Engineers, the City of Tulsa, W.M. Matthews, T. Clyde, R.M Zamor, P. Koenig, and R. West for support, assistance, and data for Lakes Texoma, Thunderbird, Grand, Eucha, and Spavinaw. JH acknowledges support from the ERDF/ESF project Biromanipulation as a tool for improving water quality of dam reservoirs (No CZ.02.1.01/0.0/0.0/16_025/0007417). BL acknowledges support from the FA-UNIMIB for long-term monitoring of Lake Iseo. EBM acknowledges support from the UK Natural Environment Research Council funding for the long-term monitoring on Blelham Tarn and the staff of the Freshwater Biological Association and UK Centre for Ecology and Hydrology for carrying out the work. AP and JAR acknowledge support from the Ontario Ministry of the Environment, Conservation and Parks for providing data from south-central Ontario lakes ("Dorset lakes") and staff and students at Ontario's Dorset Environmental Science Centre for data collection and analysis. MR acknowledges the International Commission for the Protection of Italian-Swiss Waters (CIPAIS) for funding long-term research on Lake Maggiore. EST acknowledges Prof. Lauren Chapman (McGill University) and her team for the long-term data collection in Lake Nkuruba. MS acknowledges the City of Zurich Water Supply and the cantonal agencies of the cantons of Bern (AWA, Gewässer- und Bodenschutzlabor), Zurich (AWEL), St. Gallen (AFU), and Neuchatel for providing data for the Swiss lakes, CIPEL and INRA for data from Lake Geneva, and IGKB for data from Lake Constance. WT acknowledges support from the Belgian Science Policy Office through the research project EAGLES (CD/AR/02A) on Lake Kivu. PV
acknowledges support from the councils of the regions of Waikato, West Coast, and Bay of
Plenty for long-term sampling of lakes Taupo, Brunner, and Tarawera. KY acknowledges
support from the Clark Foundation for long-term monitoring of Otsego Lake and past and current
members of SUNY Oneonta BFS for sampling. KS and JS acknowledge the National Park
Service, W Gawley, Acadia National Park for providing data for Jordan Pond, Bubble Pond, and
Eagle Lake in Maine. The views expressed in this article are those of the authors and do not
necessarily reflect views or policies of funding agencies.

**Authorship contributions:** SFJ and KCR designed the study, compiled the data, conducted
analyses, and drafted the manuscript. GJAH, BMK, PRL, JLM, RLN, RMP, JTS, CEW and RIW
helped design the study and conduct analyses, contributed data, and edited the manuscript. All
other authors contributed data, edited the manuscript, or both.

**Author information:** Derived statistics used in our analyses are publicly available via the
Environmental Data Initiative (EDI) repository at:
https://doi.org/10.6073/pasta/ac8b05bb0da19032b3df3efc21f83874.

Reprints and permissions information is available at www.nature.com/reprints. The authors
declare no competing interests.

Correspondence and requests for materials should be addressed to KCR (rosek4@rpi.edu).

377
379 **Fig. 1 | Trends in dissolved oxygen and temperature.** a-c, Density plots of trend magnitudes
380 for a temperature (°C decade⁻¹), b DO concentration (mg L⁻¹ decade⁻¹) and c DO percent
381 saturation (% decade⁻¹). Red distribution indicates surface water trends and blue indicates deep-
382 water trends. The x-axis range for each plot covers two standard deviations from the median, or
383 approximately 95% of data. Vertical dashed lines indicate median trends, and the zero trend is
384 highlighted with a black vertical line.

385
386 **Fig. 2 | Solubility effects and changes in temperature and DO concentration through time.**
387 **a, b,** Observed vs. predicted change in DO concentration (mg L⁻¹) due to solubility for surface
388 (**a**) and deep (**b**) waters. Solid black line is the 1:1 line and the blue line is loess smoothed, while
389 the gray regions are 95% confidence intervals. **c-f,** Smoothed curves of GAMM models, showing
390 deviation from the mean model predictions for selected response variables with year as the
391 predictor variable. Gray regions represent one standard error from the predicted line for **c,**
392 temperature (°C) and **d,** DO (mg L⁻¹) through time for surface waters. The red line represents
393 lakes where both surface temperature and DO were increasing (n = 87) and the blue line is all
394 other lakes (n = 332). **e,** Temperature and **f,** DO for deep waters.

395
396 **Fig. 3 | Interaction of productivity and temperature in surface waters.** a, Predicted
397 probability of a lake having both increasing surface temperature and DO concentration from a
398 fitted logistic regression model at three different mean surface water temperatures: 21°C (blue),
399 25°C (black), 28°C (red) b, Predictions of DO trends from a fitted multiple regression model for
400 chlorophyll (used as a surrogate for primary productivity) at these same temperatures (legend
401 same as a) c, The interaction of water clarity (measured as Secchi depth in m) and surface-water
402 temperature (°C) and their effects on surface DO (mg L⁻¹) from fitted generalized additive mixed
403 models (GAMM) d, Most lakes exhibited increasing surface temperatures and decreasing DO
404 concentration consistent with solubility effects, but a subset of lakes (n = 87) have both
405 increasing surface temperature and DO concentration.

406

407 **Fig. 4 | Effect of changes in water clarity and density difference on deep-water DO saturation change.** a, Change in % saturation versus change in water clarity (Secchi depth). b, Change in % saturation versus change in water column density difference between surface and deep waters. The number of lakes in each quadrant in a and b are indicated by text. c, Predictions of change in % saturation from a fitted multiple regression model for change in water clarity at three density changes. d, Predictions of change in % saturation from a fitted multiple regression model for change in density difference at three clarity changes. Note that for both c and d the origin sits at no change in either predictor.

415
<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /><figcaption>b</figcaption></figure>

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

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

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

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

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>b</figcaption></figure>

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

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

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

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

Change in
density
difference
(kg m⁻³)

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

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

Change in
Secchi
depth
(m)
416 **Supplemental information**

417 There are seven supplemental information tables and four supplemental information
418 figures. Tables S1 and S2 are referenced in text. Table S3 describes data contributors for this
419 project and Table S4 provides location and trend information for each lake. Trend data were not
420 reported for a) two lakes where providers did not provide permission to publish data but that
421 were included in trend analyses (Annecy and Geneva; ‘NP’ in table S4), b) lakes had less than 15
422 years of data at a given depth (not shown in table), or c) deep-water trends in lakes that did not
423 thermally stratify (‘NA’ in table S4). In one lake (T Bird), epilimnetic water was artificially
424 aerated and this depth layer was excluded from analysis. Table S5 presents statistics associated
425 with spatial autocorrelation analyses. Table S6 describes trends over the entire population of
426 lakes versus a sub-sample of lakes after accounting for the large numbers of samples obtained in
427 lake-rich regions. Table S7 describes trends and uncertainty in trends over two time periods for
428 subsets of lakes having data for at least 80% of years: 1980-2017 and 1990-2017. Fig. S1
429 presents the results of GAMM analysis of trends zoomed out to visualize distribution of residuals
430 for surface and deep-water temperature and dissolved oxygen trends. Fig. S2 presents the partial
dependency plots for the top predictors of changes in deep-water DO percent saturation as
431 determined by a random forest analysis. Fig. S3 presents partial dependency plots for the top
432 predictors of changes in water column density difference between surface and deep waters as
433 determined by a random forest analysis. Fig. S4 presents the locations of lakes used in this study
434 (n=393).
435
436
437 **Figure S1** | Results of GAMM analysis of trends zoomed out to visualize distribution of

438 residuals. **a**, Surface-water temperature (°C) **b**, Deep-water temperature (°C) **c**, Surface-water

439 DO (mg L⁻¹) and **d**, Deep-water DO concentration (mg L⁻¹).

440

441 **Figure S2** | **a-f**, Partial dependency plots from a random forest algorithm of deep-water change

442 in % dissolved oxygen saturation (Δ Sat) in the last five years of record relative to the first five

443 years of record for each lake. Plots are ordered by predictor variable importance, decreasing in

444 importance from the upper left to lower right (a to f). Vertical red lines indicate zero change in

445 predictor variable and hash marks on the x-axis indicate lake distribution deciles. Partial

446 dependencies indicate the relationship between predictor and response variables when holding

447 other variables at their mean value. Lakes that experienced no change in either water clarity or

448 density difference between surface and deep waters exhibited little change in deep-water

449 saturation (see also, Fig. 4).

450

451 **Figure S3** | Drivers of the change in density difference between surface and deep waters. **a-f**,

452 Partial dependency plots from a random forest algorithm of deep-water change in water column

453 density difference in the last five years of record relative to the first five years of record for each

454 lake. Plots are ordered by predictor variable importance, decreasing in importance from the

455 upper left to lower right (a to f). Vertical red lines indicate zero values for predictor variable and

456 hash marks on the x-axis indicate lake distribution deciles. Partial dependencies indicate the

457 relationship between predictor and response variables when holding other variables at their mean

458 value.

459
460 **Fig. S4** | Locations of lakes used in this study (n=393).

461
Surface Temperature

a

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

Deep Temperature

b

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

Surface DO

c

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

Deep DO

d

<figure><img src="image_19.png" /></figure>
Change in percent saturation (Δ Sat)

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

Change in Secchi Depth (m)

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

Change in density difference (kg m⁻³)

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

In (deep layer thickness (m))

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

Deep water DO saturation (%)

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

Deep DO concentration (mg L⁻¹)

<figure><img src="image_25.png" /><figcaption>f</figcaption></figure>

Change in deep temperature (°C)
<caption>Change in density difference (kg m⁻³)</caption>

a

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

Change in deep water temperature (°C)

b

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

Summer wind speed trend (m second⁻¹ decade⁻¹)

c

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

Summer air temperature trend (°C decade⁻¹)

d

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

Mean winter air temperature (°C)

e

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

Winter air temperature trend (°C decade⁻¹)

f

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

Shortwave radiation trend (W m⁻² decade⁻¹)
<figure><img src="image_32.png" /></figure>
462 **Methods:**

463 **Overview**

464 Our methods here describe how we 1) compiled and quality-checked data, 2) interpolated
465 and delineated water layer strata, and 3) statistically analyzed these data. Our statistical analyses
466 focused on characterizing long-term trends in climate characteristics (air temperature, wind
467 speed, precipitation, and short-wave radiation), DO concentration and saturation, water
468 temperature, and deep-water habitat quality; identifying and characterizing potential non-
469 linearity in DO concentration and water temperature through time; characterizing the relationship
470 between DO concentration changes and solubility, chlorophyll, and land use; identifying the
471 predictors of changes in deep-water DO saturation, and characterizing meteorological drivers of
472 surface temperature trends. These methods are described in detail in the sections below.

473 **Data compilation and quality control**

474 We compiled lake temperature and DO concentration water column measurements from a
475 wide range of government, university, and not-for-profit sources (Fig. S4 and Tables S3 and S4).
476 To assess long-term trends in temperature and DO concentration, we required profiles be made at
477 least once annually during the peak summertime stratification (defined as the late summer
478 period, July 15 - August 31 for northern hemisphere lakes and January 15 - February 28 for
479 southern hemisphere lakes) offshore (e.g., nearest the deepest location in each lake) for at least
480 15 years. In some larger lakes (n = 6 lakes), we used profiles from two separate locations if the
481 lake had more than one distinct basin and treated these as separate waterbodies. For some
482 analyses other than long-term trend analyses we included lake time series data less than 15 years
483 long, but always at least 10 years in duration (described below).
We conducted quality control on the compiled data as follows. We first removed impossible values, defined as those outside the range 0-40 for both temperature (units: °C) and DO concentration (units: mg L⁻¹). We then removed profiles from consideration if our initial quality control step process removed greater than 95% of the profile or if the profile had less than three distinct depth points. To reduce the potential impacts of DO measurements made when sensors sat on or in sediments, we removed the deepest measurement for individual profiles if the maximum depth for that profile exceeded the maximum depth of 90% of the remaining profiles for a given lake.

Not all profiles surveyed the entire water column. Some lakes had some profiles where the shallowest depth was greater than 0 (meaning near-surface measurements were not made), yet temperature measurements showed the nearest surface measurements were within the epilimnion. In these cases, we made the assumption of uniform DO and temperature from the surface to the shallowest measurement and added a 0 m depth point. We did this by either 1) changing the minimum depth in the profile to 0 if it was less than 0.5 m, 2) adding a 0 depth point and assigning temperature and DO values equal to that of the minimum depth point if the minimum depth point was greater than or equal to 0.5 m but less than or equal to 3 m. If the minimum depth was greater than 3 m, we excluded the profile from analyses. If there were multiple values of either temperature or DO for a given depth, the mean value at this depth was used. These operations and all further analyses were conducted in R version 3.4.2²⁹.

In total, the above QA steps removed 2,040 profiles out of a total of 25,023 (8.2%). After processing and removing eight non-temperate lakes, we had 22,574 DO profiles with corresponding temperature profiles. There was a median of 2.1 profiles per year (range: 1-38) and 23 years of data per lake (see also, Table S4).
507 Profile interpolation and strata delineation

508 In order to generate a dataset with consistent depth resolution within and among lakes,
509 we interpolated each temperature and DO profile from depth 0 m to the deepest depth of each
510 profile at intervals of 0.5 m using the pchip function of the R package pracma³⁰. This
511 interpolation procedure preserves the overall shape of the profile by preventing overshooting of
512 data values³⁰. Following interpolation, we calculated temperature and stability characteristics
513 using the R package rLakeAnalyzer³¹. We delineated the epilimnion and hypolimnion using the
514 meta depths function (slope = 0.1, seasonal = FALSE), which calculates the top and bottom
515 depths of the metalimnion³¹. If the range of temperatures through the profile is less than 1°C, the
516 meta depths function does not return values for the metalimnion (i.e., the profile is not
517 considered stratified).

518 Many lakes did not have a well-defined hypolimnion. To identify those with a
519 hypolimnion, we first removed lakes where the meta depths function failed to calculate a bottom
520 metalimnion depth for more than 10% of profiles. We then calculated the mean of the maximum
521 profile depths across all profiles for each lake, to get a mean profile depth for the lake. If the
522 mean value of the bottom of the metalimnion for a lake was shallower than the calculated mean
523 profile depth for that lake, it was considered to have a hypolimnion. We defined "surface waters"
524 as all depths shallower than or equal to the top metalimnetic depth and "deep waters" as all
525 depths deeper than the bottom depth of the metalimnion.

526 Characterizing trends in dissolved oxygen and temperature

527 We calculated the mean surface- or deep-water temperature and DO concentration and
528 percent saturation. For each lake, we calculated the mean of surface- or deep-water DO
concentration or temperature for all profiles in a given year (in our defined late-summer period) to obtain a mean annual value. We then calculated the percent DO saturation from temperature, DO concentration, and lake elevation data³². Mean annual surface- and deep-water temperature and DO concentration measurements were then used to calculate long-term trends for surface waters (n = 393 lakes; median number of years per lake: 24) and deep waters (n = 260; median number of years: 24). All trends were calculated using the nonparametric Sen's slope in the R package openair³³. For trend analysis, we only used lakes with at least 15 years of data.

For deep-water trends, lakes that were essentially anoxic (average hypolimnetic DO < 0.5 mg L⁻¹) had trend magnitudes that clustered near 0 relative to other lakes. This was not unexpected as lakes with essentially no hypolimnetic DO have little potential to lose additional DO. When calculating median trends and for graphical depiction of trends (Fig. 1), we removed these lakes (n = 69; difference = 191).

We conducted several analyses to examine the potential of variability in lake data through time (i.e., not all lakes sampled all years of observation) or variability in space (i.e., some regions sampled much more heavily than others) to influence overall population level trends (see following sections and Tables S5-S6).

## Spatial autocorrelation and effects of lake clusters

Because the lakes included in this study were not uniformly dispersed over all continental land masses, we examined the potential of large numbers of lakes in relatively concentrated regions to drive overall patterns. To do this, we first examined spatial autocorrelation in trends in lake temperature and dissolved oxygen concentration using Moran's I in the R package lctools³⁴, ³⁵. This statistic ranges from -1 for data that are perfectly dispersed to +1 for data that are
perfectly autocorrelated. Values near zero suggest randomly distributed data. We observed weak
but significant spatial autocorrelation in some variables (Table S5; Moran's I values ranging 0.02
to 0.27).

Following this analysis, we examined the potential for the large numbers of lakes in some
regions to dominate overall trends we reported. We tested for potential bias by examining trends
for a subset of lakes. We identified four regions in the US with high numbers of lakes (Maine =
113 lakes, New Hampshire = 38 lakes, Missouri = 41 lakes, and Minnesota = 84 lakes). For each
of these clustered regions, we randomly subsampled 10% of the lakes. After this random
subsetting, we found that the overall trends are similar to the trends obtained from all lakes (see
Table S6). These results demonstrate that our observed population-level trends are not driven
solely by trends observed in our lake-rich regions. While our analysis focuses on temperate
lakes, we obtained data from a small number of non-temperate lakes (n=8). Including these non-
temperate lakes in our analysis (Table S6) did not change our overall results.

## Uncertainty estimates and temporal variation in trends

We conducted an analysis to compare trends, confidence intervals, and significance of
trends over two time periods: 1980-2017 (n = 80) and 1990-2017 (n = 197) to assess whether
different lake observation years influenced the overall trends in DO concentration and
temperature we observed. For each time period, we used a subset of lakes that had data for at
least 80% of years within the defined time period. Following established methods¹⁸, we
calculated a yearly anomaly in temperature and dissolved oxygen for each lake as the difference
between each year's observation and the long-term mean. We then averaged these anomalies
across all lakes and used linear regression to calculate the slope, significance, and confidence
intervals of these averaged anomalies (Table S7).
574 **Characterizing trends in climate characteristics**

575 We examined trends in air temperature, total precipitation, wind speed, and shortwave
576 radiation using the ERA-5 reanalysis from the European Centre for Medium-Range Weather
577 Forecasts (ECMWF)³⁶. This data set provides a single gridded global product with a resolution
578 of 0.25° latitude by 0.25° longitude over the period 1979-2019 available as monthly averages (air
579 temperature, wind speed, and shortwave radiation) or totals (precipitation). We used ECMWF
580 time-series data from the gridded location closest to each lake and over the two-month period
581 around when lakes were sampled (July-August for Northern hemisphere lakes, January-February
582 for Southern hemisphere lakes). We calculated temporal trends in mean summer air temperature,
583 mean summer wind speed, summer total precipitation, mean summer shortwave radiation, mean
584 winter air temperature, mean spring air temperature, mean fall air temperature using the same
585 methods we used to examine lake temperature and DO trends (see above). We then conducted a
586 multiple regression analysis to assess which of these predictor variables (trends in air
587 temperature, total precipitation, wind speed, or shortwave radiation) best explained surface-water
588 temperature trends.

589 **Trends in climatic variables over the temperate zone**

590 We delineated gridded latitude and longitudes at 2° intervals across the entire temperate
591 zone over land masses only as well as over large regions, including Asia (defined by longitude ≥
592 29.3°; latitude 23.5° to 60°) Europe and North America (longitude < 29.3°; latitude 23.5° to 60°),
593 South America and western Africa (longitude < 0°; latitude ≤ -23.5° to -60°); and southern
594 Africa, Australia, and Oceania (longitude ≥ 0°; latitude -23.5° to -60°). We then used data from
595 the ERA-5 reanalysis (see ‘Characterizing trends in climate characteristics’ in Methods for
596 details) to calculate trends in climate variables over each of these regions (Table S1).
597

Multiple regression analysis of drivers of surface water temperature trends

598

599

600

601

602

603

604

605

606

607

We conducted a multiple regression analysis of the meteorological drivers of observed surface water temperature trends. Predictors in the analysis included: summer air temperature trend, summer total precipitation trend, summer wind speed trend, summer shortwave radiation trend, winter air temperature trend, spring air temperature trend, fall air temperature trend, and mean winter temperature (as a proxy for ice cover¹⁸). We z-score standardized all variables to facilitate comparison of model coefficients across variables having different units³⁷. We verified that multicollinearity was not a problem by checking that the variance inflation factor was well below ten for all variables³⁸. We used the leaps R package to select subset models including all predictors and two-way interactions, and selected the fitted model having the lowest AIC³⁹.
Coefficients and p-values for the selected model appear in Table S2.

608

Characterizing trends in deep-water habitat quality

609

610

611

612

613

614

615

616

617

We used T<sub>DO3</sub><sup>11</sup> to quantify trends in oxythermal habitat relevant for cold-water organisms. T<sub>DO3</sub> represents the minimum temperature in the water column where DO concentration was greater than or equal to 3 mg L<sup>-1</sup> and has been used to describe habitat availability for cold-water fisheries<sup>11</sup>. To calculate trends in T<sub>DO3</sub> we excluded lakes where the DO concentration was higher than 3 mg L<sup>-1</sup> across all depths in all profiles. For the remaining lakes, we calculated T<sub>DO3</sub> for each profile. If a given profile did not have DO below 3 mg L<sup>-1</sup>, we assigned it the minimum temperature in the profile. We then calculated an annual mean T<sub>DO3</sub> for the late summer period and excluded lakes that had ≤ 15 years of data. This left 369 lakes where DO went below 3 mg L<sup>-1</sup> at least once.

618

Non-linearity in DO and temperature through time
We conducted a generalized additive mixed model (GAMM) analysis to characterize overall response of lake temperature and DO concentration through time and to identify any non-linearity. GAMMs fit a smooth function of the predictor variables showing the relationship of the predictors to the response variable⁴⁰. We conducted separate analyses for four response variables, surface-water temperature, surface-water DO concentration, deep-water temperature, and deep-water DO concentration. For each GAMM, our only predictor variable was the year, resulting in models that show the change in the response variable through time. We used the gamm4 function of the gamm4 package to fit these models using the default thin plate spline for smooth terms⁴¹. Gamm4 uses penalized regression splines of moderate rank for the smooth function. For two of these models we used a normal error distribution. Because residuals for the deep-water temperature analysis were skewed, we used a gamma distribution. Residuals in the deep-water DO analysis were also skewed, but because there were a large number of 0 values we used a Tweedie distribution instead of a gamma distribution. We limited this analysis to data from 1970 and later and included all lakes with data in the specified time period (total lake n = 419). To account for the non-independent nature of the repeated measurements through time within each individual lake, the slope and intercept were allowed to vary randomly by lake⁴².

We next conducted a GAMM to understand how surface water DO concentration responded to temperature and productivity (n = 419 lakes). We used Secchi disk depth as a surrogate for productivity¹⁹. We included fixed effects of mean summer surface water temperature, mean Secchi depth, and the interaction of these two terms in the model. We included a random intercept and slope by year within each lake and included a corresponding year fixed effect.

## Relationship between dissolved oxygen concentration changes and solubility
To determine the relative importance of solubility in explaining changes in DO concentration, we calculated the expected change in DO concentration due to solubility alone and compared this amount to the observed DO change. To do this, we first calculated the difference between the observed mean DO concentration across the last five years and the first five years of record for each lake, requiring a minimum of ten years of data per lake (n = 415 lakes for surface (Fig. 2a); n = 259 lakes for deep (Fig 2b)). We then calculated the expected change due solely to solubility and compared observed to expected DO changes. Specifically, we calculated the mean percent saturation in the first five years by first calculating the mean DO saturation for each water column layer (surface or deep waters) and then calculated the mean of all of these values. We then used an analogous approach to calculate mean temperature, DO concentration, and mean DO concentration at 100% saturation in the last five years of record for each lake. Once we calculated these values, we multiplied the mean DO concentration at 100% saturation by the decimal value of percent saturation in the first five years of record. This product represents the expected DO concentration if the percent saturation in the last five years of record remained the same as it was in the first five years of record. In other words, we removed the effect of temperature so that if all changes were due solely to solubility, observed changes in DO concentration would be identical to this value.

## Relationship between dissolved oxygen trends and chlorophyll

We used multiple regression to test if chlorophyll concentration and surface-water temperature were predictors of lakes having both increasing surface DO concentration and temperature trends. We first calculated the long-term mean late-summer surface-water (epilimnetic) chlorophyll concentration, which was available for 162 lakes having at least ten years of chlorophyll measurements. We next predicted DO concentration trends using
chlorophyll and mean surface-water temperature as independent variables. We first fit the linear regression models, starting with a full model that included the interaction of chlorophyll and temperature. We then fit all subset models and selected the model with the lowest AIC value⁴³. Using this selected model, we predicted DO concentration trends at three different mean epilimnetic temperatures (21, 25, and 28°C) across the observed values for chlorophyll.

## Relationship between dissolved oxygen trends and land use

We used logistic regression to better understand the drivers of increasing DO concentration in lakes with increasing surface-water temperatures, using land use/land cover data to model the probability of this phenomenon⁴⁴. Logistic regression predicts the probability of a binary response outcome for different values of predictor variables. Predictors in our logistic regression included the percent of agriculture and developed land cover in the watershed and the mean surface-water temperature over the last ten years of record because these land use characteristics have been associated with increased growth of some phytoplankton taxa in warmer lakes⁵, ²¹. Our binary response was: either a lake had both increasing surface temperature and DO concentration (1) or it did not (0). We tested for all two-way interactions and all main effects. We used the National Land Cover Database 2011 to derive land cover metrics for US lakes⁴⁵. We considered any land falling into any of the developed classes as developed (Developed – Open Space, Developed – Low Intensity, Developed – Medium Intensity, Developed – High Intensity). We tested the goodness of fit of the final model using the Hosmer-Lemeshow test, available in the ResourceSelection R package (function hoslem.test)⁴⁶. This test showed an acceptable goodness of fit (P = 0.166). The final number of lakes for analysis that had both land cover data and sufficient data to calculate trends was 326.

## Identifying the predictors of changes in deep-water DO saturation
We first used a random forest algorithm to obtain predictors of the observed change in percent saturation (i.e., drivers beyond pure solubility effects) in deep waters⁴⁷. We used the percent increase in mean squared error as a measure of predictor variable importance. We conducted the random forest algorithm analysis using the randomForest package⁴⁸. For each analysis, we only used lakes that had no missing values for any of the predictor variables (n = 224 lakes).

For the random forest algorithm, the response variable was the change in mean DO percent saturation in the last five years of record relative to the first five years of record for each lake (Δ Sat). A positive Δ Sat indicated an increase in percent saturation while a negative Δ Sat indicated a decrease in percent saturation. Predictor variables included mean hypolimnetic DO percent saturation, DO concentration, temperature, and thickness of the hypolimnion (ln transformed), mean Secchi depth, ln of mean lake depth, log10 of residence time, change in hypolimnetic thickness, change in hypolimnetic temperature, change in Secchi depth, and change in the density difference between surface and deep waters. Mean lake depth and residence time were obtained from the HydroLakes Database⁴⁹. We calculated the density difference across the water column using rLakeAnalyzer to calculate densities for each interpolated depth point in each water column profile³¹. If a given profile was stratified, we then used the mean epilimnetic density and the mean hypolimnetic density and calculated the difference between these densities. If a given profile was not stratified, we took the mean density across the top two meters and the mean density across the bottom two meters and calculated the difference between these densities. We also included trends in the following ERA-5 meteorological variables: summer, fall, and winter air temperature, summer shortwave radiation, and summer wind speed. Finally we included mean winter air temperature as a proxy for ice cover¹⁸.
Following the above analysis, change in the density difference between surface and deep waters came out as an important predictor. Although this could be explained by increased surface water temperatures driven by meteorological variables, it is possible that other changes, such as water clarity²⁵, could also explain changes in density difference. To disentangle the drivers of changes in water column density differences, we conducted another RF using the same predictor variables as the above analysis but changing the response variable to the change in the density difference. We did not include the response variable from the first analysis (Δ Sat). The six most important variables are presented in Fig. S3.

Based on results of the RF analysis, we conducted a multiple regression analysis to predict change in percent saturation (Δ Sat) for different levels of predictor variables (ln of mean lake depth, change in the density difference across the water column, and change in Secchi depth). We used a subset of lakes where mean deep-water DO concentration exceeded 0.5 mg/L to avoid lakes with little potential to lose DO. Predictor variables were selected because they were the three most important variables identified by RF, except we substituted ln mean lake depth for ln deep layer thickness. This substitution was made because models using ln of deep layer thickness demonstrated substantial non-linearity in plots of residuals against fitted values. Models built with ln mean lake depth greatly improved these patterns and these two variables were correlated (r = 0.51). We first fit the multiple regression models starting with a full model that included all predictors and two-way interaction terms. We then fit all subset models and selected the model with the lowest AIC value⁴³. Using this selected model, we predicted Δ Sat at three different values of each of the two predictors change in Secchi depth (P < 0.001) and change in water column density difference (P < 0.001), with ln mean lake depth held at the median value.
734

735 **Data Availability:**

736 Many of the datasets analyzed during this study are publicly available on-line and associated
737 links can be found in supplementary Table S3. Derived statistics are publicly available via the
738 Environmental Data Initiative (EDI) repository at

739 https://doi.org/10.6073/pasta/ac8b05bb0da19032b3df3efc21f83874. Most lakes are included
740 here, but we note that due to the collaborative nature of this project and a wide range of data
741 provenance, it was not possible to include every lake in this repository. Data not otherwise
742 already publicly available are available upon request from the corresponding author pending
743 permission from the appropriate data provider.

744

745 **References:**

746 29. R Core Team. R: a language and environment for statistical computing. R foundation for
747 statistical computing, Vienna, Austria (2017).

748 30. Borchers, H. W. pracma: Practical Numerical Math Functions. R package version 2.1.5
749 https://CRAN.R-project.org/package=pracma (2018).

750 31. Winslow, L. A., et al. rLakeAnalyzer: Lake Physics Tools. R package version 1.11.4.
751 https://CRAN.R-project.org/package=rLakeAnalyzer (2017).

752 32. Winslow, L. A., et al. LakeMetabolizer: An R package for estimating lake metabolism from
753 free-water oxygen using diverse statistical models. *Inland Waters*, **6**, 622-636 (2016).

754 33. Carslaw, D. C., & Ropkins, K. Openair – An R package for air quality data analysis. *Environ. Model. Softw.*, **27-28**, 52-61 (2012).
755
34. Moran, P. A. P. The interpretation of statistical maps. *J. Roy. Stat. Soc. B Met.*, **10**, 243-251 (1948).

35. Kalogirou, S. *lctools: Local Correlation, Spatial Inequalities, Geographically Weighted Regression and Other Tools*. R package version 0.2-7. https://CRAN.R-project.org/package=lctools (2019).

36. Copernicus Climate Change Service (C3S). ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate. Copernicus Climate Change Service Climate Data Store (CDS), Accessed 10/1/2019. https://cds.climate.copernicus.eu/cdsapp#!/home

37. Gelman, G., & Hill, J. *Data Analysis Using Regression and Multilevel/Hierarchical Models*. Cambridge University Press, New York (2007).

38. Quinn, G. P., & Keough, M. J. *Experimental Design and Data Analysis for Biologists*. Cambridge University Press, Cambridge, U. K. (2002).

39. Lumley, T. *leaps: Regression Subset Selection (based on Fortran code by Alan Miller)*. R package version 3.1. https://CRAN.R-project.org/package=leaps (2020).

40. Wood, S. N. *Generalized Additive Models: An Introduction with R* (2nd edition). CRC Press. Boca Raton, FL (2017).

41. Wood, S., & Scheipl, F. *gamm4: Generalized Additive Mixed Models using 'mgcv' and 'lme4'*. R package version 0.2-5. https://CRAN.R-project.org/package=gamm4 (2017).

42. Pinheiro, J. C., & Bates, D. M. *Mixed Effects Models in S and S-Plus*. Springer, New York (2000).
776 43. Burnham, K. P., Anderson, D. R., & Huyvaert, K. P. AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. *Behav. Ecol. Sociobiol.*, **65**, 23-35 (2011).

779 44. Hosmer, D. W., & Lemeshow, S. *Applied Logistic Regression* (2nd edition). John Wiley and Sons, Inc., New York (2000).

781 45. Homer, C. G., et al. Completion of the 2011 National Land Cover Database for the conterminous United States – Representing a decade of land cover change information. *Photogramm. Eng. Remote Sensing*, **81**, 345-354 (2015).

784 46. Lele, S. R., Keim, J. L., & Solymos, P. ResourceSelection: Resource Selection (Probability) Functions for Use-Availability Data. R package version 0.3-2. https://CRAN.R-project.org/package=ResourceSelection (2017).

787 47. Cutler, D. R., et al. Random forests for classification in ecology. *Ecology*, **88**, 2783-2792, (2007).

789 48. Liaw, A., & Wiener, M. Classification and regression by randomForest. *R News*, **2**, 18-22 (2002).

791 49. Messager, 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, doi:10.1038/ncomms1360 (2016).

794