Main

The availability of dissolved oxygen (DO) is fundamental to the health and functioning of lake ecosystems but an alarming number of lakes experience rapid deoxygenation today1,2,3. Oxygen is an essential driver of the biogeochemical processes that regulate nutrient cycling, drinking water quality, biodiversity and aquatic food webs4. Aerobic respiration can deplete oxygen concentrations at depth when seasonal stratification prevents aeration of the water column5. This detrimental condition, termed hypoxia (DO < 2 mg l−1), diminishes the cold-water habitat of aquatic life, including commercially important fish6 and can cause mass mortality if persistent7. Anoxia (DO < 0.5 mg l−1) promotes the reductive dissolution of mineral-bound carbon and phosphorus8,9 which fuels harmful algal blooms10 and limits microbial oxidation of methane, a potent greenhouse gas11. The negative environmental consequences of deoxygenation ultimately compromise human livelihoods and health12.

Historically, oxygen depletion in lakes has been associated with excessive inputs of nutrients (eutrophication) from urban and agricultural sources, which increase oxygen demand via enhanced primary production and decomposition of sinking biomass1,13. Elevated loading of terrestrial organic matter (browning) linked to soil recovery from acidification, climate-driven shifts in precipitation patterns and land use change14,15,16, could similarly contribute to deoxygenation. Climate warming has the potential to cause further oxygen loss via prolonged summer stratification17,18 and enhanced aerobic respiration rates19. Seasonal hypoxia tends to develop more rapidly in stratified lakes with high nutrient concentrations20 and a large sediment surface area relative to the overlying water volume (Ased/Vwater), as benthic respiration is often the dominant oxygen sink5,21. This implies that oxygen levels in specific lake classes, such as small, eutrophic or humic lakes, may be more sensitive to climate forcing but whether predictors of seasonal hypoxia extend to long-term deoxygenation is unclear.

Most of the world’s lakes are located in boreal and arctic regions22, where the climate is warming at two to four times the global rate23. Rising air temperatures, atmospheric stilling and earlier ice loss have prolonged summer stratification in lakes across the region24,25,26 and, combined with widespread browning15,27, put pristine freshwater ecosystems at increasing risk of hypoxia28. At the same time, declining phosphorus and nitrogen stocks in northern landscapes associated with reduced terrestrial loading and atmospheric deposition29,30,31,32 could limit fresh organic substrate fuelling oxygen demand. Moreover, overall shortening of the ice-cover season33 may prolong spring mixing and enhance aeration before summer stratification34,35. The rapid accumulation of counteractive anthropogenic pressures means that the sign and magnitude of long-term bottom-water DO trends in northern lakes remains unknown but a continental-scale assessment, as was successfully done in 429 temperate lakes3,18, is lacking. Our inability to resolve the global extent and causes of deoxygenation hampers international commitments to address this dimension of climate impacts and protect freshwater quality.

In this study, we conducted a large-scale systematic analysis of long-term water quality monitoring data from northern lakes, comprising 2.6 million observations from 8,288 lakes in Fennoscandia and North America from 1960 to 202236. A map of the sampled lakes is shown in Extended Data Fig. 1. We compiled, harmonized and quality-checked measurements of DO, temperature, total phosphorus (TP), total nitrogen (TN), colour (a measure of chromophoric dissolved organic matter, CDOM) and lake morphometry from governmental databases and non-profit organizations via data portals and data requests. The goals of our analyses were to: (1) identify long-term changes in the bottom-water oxygen content of northern lakes, (2) assess whether standard predictors of seasonal hypoxia (surface area, maximum depth and trophic state) apply to long-term DO trends and (3) evaluate potential mechanisms driving long-term (de)oxygenation, including persistent shifts in stratification phenology, oxygen demand and spring aeration.

Results

Between 1960 and 2022, ice-free mean oxygen concentrations in lake bottom waters (≥0.7 × maximum depth) have decreased substantially across monitored northern lakes (Fig. 1). In ‘trend lakes’ with ≥15 year time series, the occurrence of bottom-water anoxia has increased from 39% to 61% between the first and last 5 years of measurement (n = 912). Oxygen trends scaled inversely with lake surface area from −0.40 to +0.07 mg l−1 decade−1 in lakes sized ≤10 ha and >104 ha, respectively (Fig. 2a). This scaling pattern was robust across continents (Fennoscandia and North America) and biomes (tundra, boreal forests, north-temperate forests and grasslands)37 (Extended Data Fig. 2). Because small lakes and ponds are disproportionally abundant22, our results imply that most northern lakes have experienced major deoxygenation. In our dataset, small lakes are under-represented compared to their real relative abundance. Adjusting for this bias via weighted resampling, we estimated that bottom-water DO has decreased below the critical hypoxic threshold (2 mg l−1) in over half of northern lakes (Fig. 1a, black line).

Fig. 1: Decadal time series of DO.
figure 1

a–d, Bottom-water oxygen concentrations classified by lake surface area (a), maximum depth (b), trophic state (c) and humic content (d) across 7,394–8,288 lakes, depending on data availability. Annual class-medians (dots) were computed from ice-free season lake–year mean DO corrected for bias toward mid-season sampling (Methods). GAM smooths were weighted by inverse variance to account for uncertainty at low sample counts. Shaded areas represent 95% confidence intervals (CIs). The black curve in a shows a weighted population-median estimate based on the size distribution of Swedish lakes (Methods). In a and b, legend labels define class boundaries (for example in a, green marks Asurf = 102–103 ha). Trophic state categories are oligotrophic (O), mesotrophic (M) and eutrophic (E). Humic content categories are clear (C), humic (H) and polyhumic (PH). Dashed and dotted lines mark hypoxia (DO < 2 mg l−1) and anoxia (DO < 0.5 mg l−1), respectively.

Full size image
Fig. 2: Trends in DO and stratification strength.
figure 2

a–c, Box plots of long-term trends (≥15 years) of bottom-water DO (a), stratification strength (Δρ) (b) and sediment surface area to water volume ratio (Ased/Vwater) (c) by logarithmic lake surface area class. d, GAM residuals of DO trends as a function of Ased/Vwater, Δρ trends and sampling day-of-year trends. Boxes bound the interquartile range (IQR) (25th to 75th percentile) and whiskers extend to the smaller quantity of data extremes or medians ± 1.5× IQR. Dots and lines mark means and medians, respectively. In a, b and d, open symbols signify means whose 95% CIs include 0. Outliers are shown by grey dots. Distinct letters (top) mark statistically different groups (Padj < 0.05, Conover–Iman tests). Bottom labels show the sample size (number of trend lakes) in each surface area class, which differed per variable depending on data availability.

Full size image

The areal scaling of long-term bottom-water oxygen trends (Fig. 1a) divided into four distinct size groups (Conover–Iman test, Padj < 0.05) at 102, 103 and 104 ha (trend lakes, Fig. 2a), with significant downward trends in lakes <103 ha and non-significant mean increases in the largest lakes (>104 ha). Among lakes classified by their maximum depth (zmax, 20th percentiles), those of intermediate zmax (4–31 m) showed substantial oxygen declines, whereas bottom waters of shallow (zmax < 4 m) and deep (zmax ≥ 31 m) lakes showed no clear long-term trend and generally remained oxygenated throughout the ice-free season (Fig. 1b). In addition, long-term oxygen loss was markedly more pronounced in eutrophic and polyhumic (highly coloured) lakes compared to less productive and clearer water bodies (Fig. 1c,d). These covariates are consistent with empirical models of seasonal oxygen loss as a function of basin morphometry and nutrient concentrations20 and demonstrate that similar emergent mechanisms modulate long-term oxygen trends.

To identify the dominant drivers of the observed oxygen trends we conducted random forest (RF) analyses38 of ice-free mean bottom-water DO saturation (DOsat) and its long-term trend. Separating the integration timescale is necessary because cross-scale interactions among drivers identified in seasonal and comparative studies may shift their importance on multiyear scales39. Predictor variables included stratification strength (Δρ), sampling day-of-year, temperature, colour, TP, TN and their respective long-term trends and Ased/Vwater. We used surface-water colour and nutrients because in the hypolimnion, anoxia-induced internal loading confounds causal relations with DO8. The RF models explained 76% and 28% of the variance of DOsat and DOsat trends, respectively. Predictor variable importance scores ranked stratification strength as the primary covariate of seasonal and long-term DO depletion, followed by Ased/Vwater (Extended Data Fig. 3). Partial dependence plots show that bottom-water DOsat tended to decrease with increasing Δρ, Ased/Vwater and surface-water TP and TN (Fig. 3a). DOsat trends were more negative (deoxygenation) in lakes experiencing strong positive Δρ trends and surface-water browning and eutrophication (Fig. 3b). To assess changes in oxygen consumption, we computed the seasonal bottom-water oxygen depletion rate (∂DO/∂t) from annual sets with more than four observations of decreasing DO concentrations at each profiling location under stratified conditions (Δρ ≥ 0.1 kg m−3) (n = 2,528). Assuming negligible transport of DO through the thermocline40, this rate closely approximates biogenic volumetric oxygen demand5. We found that ∂DO/∂t trended negative in 64% of a small but representative set of trend lakes (−1.99 ± 1.31 µg−1 l−1 d−1 decade−1, n = 53), suggesting that DO consumption rates have not systematically increased in northern lakes and therefore do not generally explain the trend towards hypoxia.

Fig. 3: Environmental drivers of DO saturation.
figure 3

a,b, RF partial dependence plots of ice-free mean bottom-water DOsat (a) and its long-term trend (b). a shows the six most important predictor variables identified in the RF model; Δρ, Ased/Vwater, sampling day-of-year and surface-water colour, TP and TN. b includes the top three predictors of the DOsat trend (Δρ trend, Ased/Vwater, sampling day-of-year trend) as well as trends in surface-water colour, TP and TN. Bars summarize statistics of the independent variables and follow the legend of Fig. 2. Numbers rank the permutation importance z-score relative to that of the most important feature. The y axes are the same scale for each subpanel; grey dashed lines mark the intercepts. The complete importance ranking of covariates is shown in Extended Data Fig. 3.

Full size image

We evaluated long-term trends in bottom-water oxygen during the spring overturn and the start of summer stratification (first 3 weeks after ice-off). Size-class-mean concentrations varied between 5.6 and 10.6 mg l−1 (48%–84% saturation) and increased with lake size (Fig. 4). Contrary to summer observations, initial oxygen trended positive in all but the smallest lakes. To examine how these partly opposing trends are linked, we estimated decadal changes in seasonal DO cycles by fitting a hierarchical generalized additive model (HGAM)41, using a normalized time variable—day-of-year scaled to mean ice-on and ice-off dates—to account for geographical differences in the length of the ice-free season (Methods). The HGAM revealed that long-term trends toward oxygenation in spring can offset deoxygenation from earlier onset of stratification (Fig. 5) and that the balance between the two processes is, in part, a function of lake morphometry. Enhanced spring aeration elevated summer DO in some of the largest lakes with negligible shifts in stratification phenology (≥104 ha) but the net offsetting effect diminished with increasing thermal stability trends toward mid-sized lakes (Figs. 2b and 5).

Fig. 4: Decadal time series of DO concentrations in spring.
figure 4

Annual median bottom-water oxygen concentrations in the first 3 weeks after ice-off by logarithmic lake surface area class. GAM smooths were weighted by inverse variance to account for uncertainty at low sample counts. Shaded areas represent 95% CIs. Legend labels define area class boundaries (for example, green marks Asurf = 102–103 ha). Dashed and dotted lines mark hypoxia (DO < 2 mg l−1) and anoxia (DO < 0.5 mg l−1), respectively.

Full size image
Fig. 5: Shifts in the seasonal cycle of DO.
figure 5

HGAM estimates of mean bottom-water DO as a function of (scaled) day-of-year, year, lake surface area (Asurf), maximum depth (zmax) and relative sampling depth (zsample/zmax) (Methods). Predictions are shown for years 1970 and 2020 at a depth of 0.85 × zmax in order-of-magnitude lake surface area classes and class-median zmax (that is, for Asurf = 10 ha, the median zmax of lakes 5–50 ha). Day-of-year was scaled to lake-mean ice-on and ice-off derived from the ERA5-Land ‘lake ice depth’ product61 to account for geographic differences in ice-free season length. Vertical dashed lines mark local oxygen maxima and approximate the timing of vernal and autumnal lake mixing. Shaded areas cover 95% Bayesian CIs.

Full size image

We examined whether any one mechanism dominates the morphometric scaling of DO trends (Fig. 2a). First, trends in stratification strength (Fig. 2b) and duration (Fig. 5) increased with decreasing lake size. Second, the predictive power of Ased/Vwater in the RF analyses on both seasonal and interannual timescales (Fig. 3) suggests that the response to climate perturbations (for example, an additional day of stratification per year) is proportional to areal extent of the sediment oxygen sink (Ased/Vwater), which scales inversely with lake surface area (Fig. 2c). Third, DO at the start of the ice-free season trended positive in large lakes (Figs. 4 and 5). We constructed a generalized additive model (GAM) of DO trends with terms for Δρ trends (1), ln(Ased/Vwater) (2) and ice-free mean sampling day-of-year trends (3), plus an interaction term (1) × (2) for lakes where both DO and Δρ trends could be computed (n = 787 lakes). Term (3) accounted for progressively later sampling in some lakes which could have biased DO trends low (Fig. 3b) but was not significant (P = 0.112). Figure 2d presents summary statistics of the model residuals. The model explained class-mean differences between DO trends (Kruskal–Wallis test on model residuals, P = 0.19, d.f. = 4) (Fig. 2d), while significant differences remained in the residuals of GAMs excluding either term (1) (P = 0.0071, d.f. = 4) or term (2) (P < 0.0001, d.f. = 4). Spring DO trends were omitted from the GAM for lack of trend lakes (n = 84) but accounted for the non-significant, positive mean DO trend in lakes >104 ha (Fig. 2a) remaining in the model residuals (Fig. 2d). Thus, the lake-size dependency of bottom-water DO trends emerged from a synergy between biotic and physical drivers.

Seasonal anoxia is known to enhance internal loading of CDOM and P (via reductive desorption) and N (as NH4 via nitrification) in the hypolimnion of stratified lakes but the long-term effects of the observed proliferation of anoxia (Fig. 1) remain unknown. We computed seasonal rates of change of bottom-water colour, TP, TN and NH4 under distinct oxygen regimes for every lake–year with four or more observations in the stratified period (Δρ ≥ 0.1 kg m−3). As expected, we found significant accumulation of all four substances under anoxic conditions (Fig. 6). Rate magnitudes decreased with increasing DO and were negligible in highly oxic (>8 mg l−1) waters. The seasonal increase in TN under anoxic conditions was accounted for by NH4. We next computed bottom-water colour, TP, TN and NH4 trends (≥15 years) under the same oxygen regimes to assess their long-term impacts. Overall, mean bottom-water colour trends were positive, TP trends were negative and TN and NH4 trends tended to zero (Fig. 6). The magnitudes of colour trends increased with decreasing DO and were approximately four times higher in anoxic compared to oxic hypolimnia. In contrast, and perhaps counterintuitively, we found that recurring anoxia was associated with dispersed trends in bottom-water TN and NH4 and significant downward trends in TP.

Fig. 6: Seasonal rates and long-term trends of colour, N and P across oxygen regimes.
figure 6

a–d, Estimates of ice-free seasonal rates (within years, orange) and long-term trends (across years, cyan) of bottom-water colour (a), TN (b), TP (c) and NH4 (d) within distinct oxygen regimes: anoxic (≤0.5 mg l−1), hypoxic (0.5–2 mg l−1), low oxic (2–5 mg l−1), mid oxic (5–8 mg l−1), high oxic (>8 mg l−1) and ‘All’ (no oxygen criterion). Long-term surface-water trends are shown for comparison as open symbols. In the y axis labels, f stands for frequency (d−1 and yr−1). The category ‘All’ includes all bottom-water samples. Error bars mark 95% CIs of the means. Statistics associated with this figure are summarized in Extended Data Table 1.

Full size image

Discussion

Long-term deoxygenation of northern lakes was strongly coupled to increased density stratification (Δρ) (Fig. 3), which in turn is driven by climate warming18,24,25. While long-term patterns in air temperature and wind speed were near-identical across lake-size classes (Extended Data Fig. 4), smaller lakes experienced the greatest prolongation of summer stratification (HGAM, Fig. 5) and Δρ trends (Fig. 2b). This suggests that lake morphometry mediates thermal responses to climate warming. This notion is corroborated by model and field studies; in small dimictic lakes with limited fetch, reduced wind mixing confines heating to a shallow surface layer, increasing the thermal gradient42, whereas in lakes with a large surface area, wind mixing distributes heat over greater depths43 and interannual variation of stratification duration is governed by wind speed more than air temperature44.

Our results further indicate that prolonged thermal stratification, rather than fewer mixing events, is the main driver of long-term oxygen loss. Local oxygen maxima in the modelled seasonal DO cycles (Fig. 5) revealed that spring and autumn mixing occurred significantly earlier and later in the year, respectively, consistent with reports of widespread shortening of lake ice-cover periods across the northern hemisphere33. Our findings are also in agreement with previous studies in temperate lakes which report coupled increases in stratification duration and bottom-water hypoxia17,18,34,45 and, more rarely, negligible DO trends where the stratification period has not lengthened but shifted46. Here, dividing long-term DO trends by seasonal ∂DO/∂t in lakes where both estimates were available (n = 251), we estimate that an extra (median) 1.7 days of DO drawdown per decade accounted for the DO trends. This is comparable to model estimates of historic (1960–2005) prolongation of summer stratification in northern hemisphere lakes larger than 100 ha (ref. 24), about 2 d decade−1.

In addition to duration of stratification, the degree of aeration at the start of the stratified period in spring can severely impact the trajectory of oxygen loss over summer5,35. We found that initial (spring) oxygen concentrations decreased predictably with lake surface area (Figs. 4 and 5). This pattern may in part be a remnant of winter oxygen levels, which, as in summer, were elevated in larger lakes47. Additionally, during the brief spring overturn in smaller lakes, their smaller fetch can limit efficient air–water gas exchange (aeration)48. The time series of initial DO (Fig. 4) and the HGAM (Fig. 5) further showed that, in spring, DO trended positive in lakes larger than ~103 ha. Such improved bottom-water aeration could indicate that earlier ice-off has caused a lengthening of the spring mixing period35 but the drivers of this widespread and consistent pattern remain unclear. We hypothesize that in large lakes with substantial thermal inertia, excess warming in winter26 has had a greater effect on the timing of ice-off than on stratification onset.

Previous studies have linked hypolimnetic oxygen loss to cultural eutrophication of surface waters and increased availability of organic matter1,10,49. In agreement with these reports, we found that long-term bottom-water oxygen trends scaled inversely with trends in surface-water TP and colour (Fig. 3b). Here, oligotrophication of northern lakes29,30,31 may have offset some of the climate-driven deoxygenation, which would be consistent with the observed diminishing of seasonal oxygen depletion rates (∂DO/∂t). Continental-scale patterns in surface-water browning14,15,27 may have exacerbated long-term bottom-water oxygen loss via enhanced heat absorption at the water surface strengthening thermal stratification50. Browning may also have increased lacustrine degradation of terrestrially derived organic matter51, although this seems inconsistent with observed ∂DO/∂t declines. Overall, biogeochemical predictors ranked consistently lower in importance than physical quantities Δρ and Ased/Vwater (Extended Data Fig. 3). Compared to climate warming, browning and nutrient trends appear to have had a limited effect on oxygen loss in northern lakes at a monitoring timescale of ∼25 years (trend lake median). Another way to interpret this result is that biologically mediated changes in lake metabolism tend to occur at a much slower pace1.

Conversely, in an increasing number of lakes experiencing annual anoxia (Fig. 1), biogeochemical cycling of carbon and nutrients may be affected long-term. The positive CDOM trends in anoxic waters (Fig. 6) imply that anoxic liberation of buried organic matter, previously shown in sediment incubations and in situ in lakes on seasonal timescales8,52,53, is not only widespread but has also increased over the last six decades. As a result, bottom-water colour trends were significantly greater than landscape-driven browning rates measured at the surface (trend lakes, pairwise one-sided Wilcoxon signed rank test, P ≪ 0.0001, n = 678), with a factor two difference between mean trends (Fig. 6). Notably, a recent study found that dissoved organic matter (DOM) ‘deprotected’ from iron complexes is highly reactive and that up to 21% can escape coprecipitation with iron upon re-aeration54. If retained in the lake, this labile DOM pool could, together with reduced substances that also accumulate in anoxic hypolimnia (CH4, NH4, S(−II), Mn(II), Fe(II))55, deplete the limited supply of epilimnetic DO mixed down during the autumn overturn56 and reinforce anoxic conditions over longer timescales. Latent oxygen sinks may be particularly important in smaller lakes with higher Ased/Vwater (Fig. 2c) and greater build-up of sediment-derived substances55. Here, the HGAM (Fig. 3) indicates that bottom-water oxygen concentrations at ice-on have declined in lakes smaller than 100 ha, but more work is needed to identify the cause(s) of this trend.

As with CDOM, we might expect prolonged stratification (and anoxia) to elevate TP concentrations long-term through increased internal loading. While bottom-water TP trends were generally more positive than surface-water TP trends (pairwise one-sided Wilcoxon signed rank test, P = 0.003, n = 820), this pattern reversed under persistently recurring anoxia (≥15 years) (Fig. 6c). We hypothesize that the strong downward TP trend in anoxic bottom waters was caused by the gradual depletion of iron-bound P reservoirs in the thin, periodically oxygenated layer of the surface sediment. One mechanism that could drive such depletion is incomplete coprecipitation of P with Fe upon re-aeration57, if the residual P is exported out of the lake or sequestered into less redox-sensitive complexes. Prolonged anoxia evidently did not exhaust the redox-sensitive organic matter fraction in surface sediments (Fig. 6a), probably because its concentration tends to be several orders of magnitude higher than that of iron-bound P58,59.

Northern lakes, which have until recently been less affected by human development, are now at severe risk of prolonged hypoxia with consequent ecological deterioration. Recent long-term studies in north-temperate lakes show the rapid diminishing of the oxythermal habitat of commercially important fish under concurrent browning and warming60 and the potential for mass kills during heat waves7. On the basis of our results, we expect a similar, ongoing impact in boreal and arctic regions, where cold-water fisheries fulfil important cultural and economic functions. Moreover, the proliferation of anoxia could have major biogeochemical implications, including increased remobilization of terrestrially derived DOM and a decrease in the amount of terrestrial carbon that is ultimately buried in lake sediments53,54. Deoxygenation has particularly affected small lakes, globally the most abundant lake type and is a widespread problem with important socio-economic consequences that demands urgent attention.

Methods

In this section, we describe how we sourced the lake water quality and morphometry data, quality-controlled the observations, computed derived quantities, such as volumetric oxygen demand and analysed the data with established statistical methods in R v.4.2.0 (ref. 62). Variables were selected on the basis of known (proxy) relations with oxygen and availability across the different datasets. The goal of the statistical analysis was to identify environmental covariates of bottom-water oxygen concentrations and their long-term trends. This approach enabled us to disentangle the impact of climate warming on deoxygenation from that of anthropogenic changes in lake metabolism associated with eutrophication and browning. A detailed description of these methods follows.

Data compilation

In this study, we focused on seasonally ice-covered lakes in the northern regions of North America (>44° N) and Northern Europe (>55° N), where most of the world’s lakes are located22 and the climate is warming significantly faster compared to the global average23. We compiled lake water quality measurements from a wide variety of governmental and not-for-profit organizations in the United States (Alaska), Canada (British Columbia, Québec, Ontario, Saskatchewan, Alberta, Manitoba, Nova Scotia, New Brunswick, Prince Edward Island and the Northwest Territories), Finland, Sweden and Norway (Supplementary Data 1). Data were accessed directly through online data portals (7,855 lakes) or obtained via data requests (433 lakes). Lake morphometry information was similarly collected from public repositories, scientific literature and data requests (Supplementary Data 2) and matched with the water quality data through common water body identifiers or manual colocation. Lakes without maximum depth or surface area information were excluded from the study. The synthesized dataset comprises 2.6 million measurements in 8,288 lakes36.

Site selection

We included freshwater lakes and ponds in the dataset but we excluded reservoirs, non-permanent lakes and actively managed water bodies. We classified managed water bodies as those undergoing physical measures against deoxygenation, such as aeration, dredging or hypolimnetic withdrawal, as well as chemical treatments to combat eutrophication which are likely to affect oxygen consumption rates, such as the addition of aluminium sulfate. Because measures to improve water quality are often costly and disruptive they are generally well documented. We identified 197 managed lakes and ponds in our dataset with help from the data providers, through news media, reports and published research and via regional authorities and local conservation organizations. Basins of large lakes that were sampled separately were treated as distinct water bodies. Each lake was assigned a new identifier to replace originals which could overlap between datasets and to merge observations from water bodies sampled by multiple organizations, for example across borders.

Variable selection and computation

Physical and chemical variables were selected on the basis of known (proxy) relations with bottom-water oxygen concentrations and availability in the different datasets. In our data synthesis, we included the following variables: water temperature, true water colour, TP, TN, ammonia and chlorophyll. TP and TN were computed as the sum of particulate and dissolved fractions (that is, unfiltered samples). When TN measurements were unavailable, such as in some profiles of the US WQP dataset63, we computed TN from total Kjeldahl nitrogen (ammonia + organic N), nitrate and nitrite. In analyses involving long-term TN trends, we confirmed that TN was computed consistently through time. We also included sampling coordinates, sampling date, sampling depth and lake altitude, surface area, mean depth and maximum depth in the dataset.

We used standard equations to compute freshwater density64 and DO saturation with respect to atmospheric concentrations65. Strength of thermal stratification was estimated as the difference between mean bottom- and surface-water density (Δρ, kg m−3). There exists no universal definition of stratification in lakes but density thresholds are frequently applied66 and we defined stratified profiles as those with Δρ exceeding 0.1 kg m−3(ref. 24). Salinity measurements were unavailable in most lakes and were not included in the computation of freshwater density. This means that in chemically stratified naturally hypersaline lakes, which can be found in select areas of the northern Great Plains in western Canada67, our RF estimate of the importance of Δρ (Extended Data Fig. 3a) is probably conservative. Trophic state classes were assigned on the basis of surface-water chlorophyll or, if unavailable, TP concentrations as oligotrophic (chlorophyll (Chl) ≤ 2.6 µg l−1, TP ≤ 12 µg l−1), mesotrophic (Chl = 2.6–7.3 µg l−1, TP = 12–24 µg l−1) and eutrophic (Chl > 7.3 µg l−1, TP > 24 µg l−1)68. Humic state was divided into three categories on the basis of surface-water colour, measured on a platinum–cobalt (Pt) colour scale; clear (<30 mg Pt l−1), humic (30–90 mg Pt l−1) and polyhumic (>90 mg Pt l−1)69. The Ased/Vwater ratio below 0.7 × zmax was estimated using hypsographs computed with the approx.bathy function in the rLakeAnalyzer R package70.

Filtering and quality control

Quality control of the synthesis dataset was conducted as follows. First, we removed observations without a sampling date or (plausible) sampling depth (0–lake maximum depth). However, bathymetric surveys based on transects can be incomplete, so in water bodies where the sampling depth consistently (in more than five profiles) exceeded the reported zmax we adjusted zmax accordingly. We standardized units of measurement across the dataset and only kept measurements with unambiguous units, for example NH4, TN or TP in molar quantities (µM) or concentrations ‘as N’ and ‘as P’ and clearly defined fractions (‘total’ or ‘dissolved’)71. We subsequently removed water quality measurements outside a predefined plausibility range, generally excluding the top 0.1% quantile: DO (0–20 mg l−1), water temperature (0–40 °C), TP (0–3 mg l−1), TN (0–60 mg l−1), water colour (0–800 mg Pt l−1) and NH4 (0–8 mg l−1). The range filter was applied to the raw data before the computation of derived quantities (see the following sections). For all analyses, except the HGAM of seasonal DO cycles (Fig. 5), we removed profiles outside the ice-free period (see section on ‘Definition of the ice-free period’). In figures comparing class-mean values, trends or rates, we removed extreme outliers (>3 s.d. from the mean) which comprised 1.3%–2.2% of values.

Definition of depth zones

We defined the bottom water of each lake as the layer below 0.7 × maximum depth and surface water as the layer above the smaller quantity of 0.5 × maximum depth and 1 m depth. A more traditional approach defines the surface and bottom layer as the epi- and hypo-limnion, respectively, but this would exclude lakes and lake-years without a clearly defined seasonal thermocline. By including periods when lakes are fully mixed, our definition does not mask the ongoing shift from polymictic to dimictic mixing regimes in many small and shallow lakes as a result of surface warming24.

Definition of the ice-free period

We estimated ice phenology for each lake from the ‘lake ice depth’ product in the ERA5-Land reanalysis dataset61. The product is based on the ECMWF weather model coupled with the 1D hydrodynamic lake model FLake72. FLake has been shown to accurately predict ice-on and ice-off dates in north-temperate lakes with a high accuracy73. In each lake, the ice-free season was defined as the portion of the period between 1 June and 31 October when the ice thickness was 0. We also computed long-term (1960–2022) mean ice-on and ice-off day-of-year for each lake as temporal reference points for the HGAM.

Determination of long-term trends

We estimated the magnitude of the interannual trends of bottom-water oxygen concentrations and environmental variables. First, we computed ice-free, depth-averaged annual mean values for each variable and lake–year. In lakes where persistent anoxia biased the oxygen trend towards 0, rather than omitting the lakes from the trend analysis entirely, we removed each middle value of three consecutive years with ice-free mean DO values below 0.5 mg l−1. This procedure preserved initial deoxygenation preceding long-term anoxia in lakes with multidecadal DO time series. We then selected the subset of lakes with time series exceeding 15 years, comprising 460–951 trend lakes depending on the variable. Omitting anoxic years reduced the number of DO trend lakes from 912 to 805. For each lake and variable, we quantified long-term changes with the Theil–Sen slope estimator74, a non-parametric technique for trend estimation that is robust to outliers and non-normality.

Determination of the seasonal oxygen depletion rate

We estimated ice-free season (within-year) seasonal oxygen depletion rate (∂DO/∂t) as the rate of change of bottom-water DO with the Theil–Sen estimator74 for time series with a minimum of four measurements per lake–year. We selected, for each lake and year, the part of the time series with decreasing DO concentrations and excluded non-stratified profiles (Δρ < 0.1 kg m−3). In this way, we excluded periods of intermittent and autumn mixing from the computation of ∂DO/∂t. The computed seasonal oxygen depletion rate reflects a balance of oxygen sources (diffusion and photosynthesis) and sinks (biochemical oxygen consumption). During summer stratification, transport of DO through the thermocline is typically orders of magnitude lower than hypolimnetic oxygen consumption40. In addition, photosynthesis is less likely to play a significant role as 87.4% of bottom-water oxygen samples were collected from the aphotic zone (<1% of surface photosynthetically active radiation at >2.7 × Secchi depth4). We therefore interpreted ∂DO/∂t to reflect primarily the volumetric oxygen demand of aquatic organisms. Even so, we cannot completely exclude an influence of oxygen sources and, as such, oxygen consumption rates inferred from ∂DO/∂t should be considered conservative.

Comparison of seasonal rates and long-term trends across an oxygen spectrum

To examine functional relations between oxygen and biochemical variables, we computed seasonal rates and long-term trends in bottom-water colour, TP, TN and NH4 under distinct DO regimes: anoxic (≤0.5 mg l−1), hypoxic (0.5–2 mg l−1), low oxic (2–5 mg l−1), mid oxic (5–8 mg l−1) and high oxic (>8 mg l−1), plus an additional class that included all DO concentrations. We excluded non-stratified profiles (Δρ < 0.1 kg m−3) to minimize the influence of turbulent mixing. For each variable, we estimated seasonal rates by first computing mean values for each sampling date and within each available oxygen regime and then estimating a rate for each lake–year and oxygen class with more than four observations. Long-term trends were similarly computed from mean values across years and oxygen classes for each lake with more than 15 years of data. Both rates and trends were computed with the Theil–Sen estimator74.

Sampling bias correction

In water quality datasets, small lakes and ponds are typically under-represented relative to their natural abundance. Lakes with surface areas <10 ha comprise over three-quarters of global lakes75 but only 26% of our dataset. We therefore computed size-class abundance weighted medians of annual lake-mean oxygen concentrations (black line in Fig. 1). Abundance fractions were derived from the Swedish national lake registry, which lists the size class of over 100,000 water bodies (Svensk Vattenarkiv (SVAR) v.2016_3, retrieved 23 September 2022 from the Swedish Hydrological and Meteorological Institute). The census is based on a compilation of maps and is considered to be a complete record of Swedish lakes ≥1 ha (ref. 75).

At the level of individual lakes, we corrected for temporal sampling biases using GAM residuals. The number of samples per month peaks in August when seasonal oxygen minima occur most frequently, which could bias ice-free mean DO estimates low. The procedure is similar to detrending of a time series using a linear regression model. Here, we fitted an HGAM of the bottom-water DO concentration with a smooth for day-of-year (lake–year means), random intercepts for each lake and a Tweedie conditional distribution to account for the disproportionate number of zeros in the response variable. We fitted the GAM with the R package mgcv41. The day-of-year term was significant (P ≪ 0.0001). Model residuals were subsequently used in pooled time-series analyses (Fig. 1). We also detected a mean tendency towards sampling later in the year (1.6 ± 0.4 d decade−1, n = 951), which we adjusted for using GAM residuals before conducting the comparative analysis of DO and Δρ trends in Fig. 2. In the RF analyses of bottom-water oxygen saturation (trends), we included day-of-year (trends) as features.

We evaluated whether the variation in Δρ with lake morphometry (Fig. 2b) could mask a geographic climate forcing pattern. Overlapping time series of ERA5 reanalysis 2 m air temperature and 10 m wind speed across lake area and depth classes rule out a bias in climate forcing (Extended Data Fig. 4) and instead confirm that lake morphometry mediates lake thermal responses to climate warming.

Random forest analysis

We conducted RF analysis38 to identify the most important environmental covariates of ice-free mean bottom-water water DO saturation (1) and its long-term trend (2). We chose to model DOsat rather than concentration to distinguish temperature effects on solubility and biochemical conversion rates. RFs are based on an ensemble of regression trees, where each tree is trained on a bootstrap sample of data and a random subset of the independent variables splits each node. Model performance is evaluated as the error of its predictions for a subset of data not used to train the model (33% of the dataset). The importance of each variable corresponds to the decrease in model performance upon shuffling of its values (permutation). RF models were fitted with the Ranger package in R76. In RF (1) we included predictor variables that represent known drivers of DO5,20; stratification strength, surface-water TN, TP and colour, bottom-water temperature, Ased/Vwater and day-of-year. In RF (2) we added long-term trends of all variables, except Ased/Vwater. We used surface-water colour and nutrients as predictor variables because, in the hypolimnion of stratified lakes, anoxia-induced internal loading may confound causal relations with DOsat. Because RF cannot easily distinguish the importance of highly correlated variables (collinearity), we computed Kendall’s τ for each independent variable pair but found none that exceeded the recommended threshold of |τ| > 0.7 (ref. 77). We used the tuneMtryFast function to determine the optimal number of splitting variables and set up the RF with 1,500 trees to ensure stable outcomes. We created partial dependence plots with the R package ‘pdp’78 to visualize the marginal effect of each predictor variable. We re-fit model (1) replacing bottom trends with surface trends of TN, TP and colour to test hypotheses about the causal nature of the correlations (see main text).

Statistical tests

To assess statistical significance of central tendency differences across classes (Fig. 2) we conducted one-way analysis of variance with the Kruskal–Wallis test79 and post hoc Conover–Iman tests for stochastic dominance with Benjamini–Hochberg adjustment for multiple comparisons80,81. The null hypothesis (no difference between groups) was rejected at Padj < 0.05. To test the hypothesis that long-term TP, TN and colour trends are greater in deep waters than in surface waters we performed pairwise one-sided Wilcoxon signed rank tests. Non-parametric tests were chosen because conditional distributions were not necessarily normal and because they tend to be more robust against outliers.

Resolving interannual changes in the seasonal oxygen cycle

We fitted an HGAM to evaluate how the seasonal cycle of bottom-water DO varies between years and across lake morphometric gradients. GAMs are a class of generalized linear models, in which the response variable depends linearly on unknown smooth functions of the predictor variable(s)41,82. Hierarchical or mixed models contain fixed and random effects to explicitly model variance–covariance structures of non-independent, clustered data (for example, oxygen observations in lakes)83. The hierarchical model allowed us to incorporate information from the full dataset, rather than a subset of lakes with long-term, high-resolution time series. The HGAM was fitted with the mgcv package41 in R v.4.2 compiled with OpenBLAS v.0.3.15, on a high-performance computing cluster hosted by the National Academic Infrastructure for Supercomputing in Sweden (NAISS). As fixed effects we included smooth terms for day-of-year (DoY), year, relative sampling depth (zsample/zmax), log(surface area) and log(maximum depth) and their two- and three-way interactions. We modelled the effects as smooth functions that were represented in the model via cyclic (DoY) and cubic (all other variables) regression spline basis expansions of the covariates. We scaled day-of-year in each lake to long-term mean ice-on and ice-off dates from the ERA5 reanalysis dataset to account for geographical differences in ice-free season length. Including random slopes and intercepts for each lake resulted in a model that was too large to fit, so we opted for lake clusters within which we could reasonably expect highly synchronous oxygen cycles. To this end, the lakes were divided into 500 groups by standardized (zero mean and unit variance) latitude, longitude, annual mean air temperature, surface area and maximum depth via k-means clustering. The number of groups was limited by the memory capacity of the computing cluster. We modelled the conditional error distribution as a Tweedie location scale family84, in which the additional distribution parameters (power (p) and scale (φ)) are also modelled as sums of smooth functions of covariates, to account for distinct conditional error distributions, for example, in the mixing periods versus at the height of summer stratification.

Seasonal cycles were predicted for years 1970 and 2020 and a relative depth of 0.85 × zmax for four order-of-magnitude lake surface area classes with corresponding class-median maximum depths, using the ‘predict’ function in mgcv. Confidence intervals of predictions were computed from standard errors based on the Bayesian posterior covariance matrix of the parameters. We omitted more uncertain predictions for sparsely sampled data regions, such as small ponds (Asurf = 1 ha) and early sampling dates (1960s). Like most linear regression models, the HGAM provides expected means at specified values of the independent variables. Predictions therefore represent the expected mean DO concentrations of all lakes in a defined variable space (for example, 10 ha lakes in 1970) and should not be interpreted as ‘typical’ patterns in specific lake types or individual lakes. We estimated the long-term change in the duration of stratification by tracking the timing of the vernal and autumnal oxygen maxima. These maxima represent model expectations of the time when most lakes of the specified morphology reach peak DO concentrations. Mean peak locations and 95% credible intervals (CIs) were computed from predictions based on repeated draws from the posterior distribution of model parameters (n = 1,000) using the fitted_samples function in the gratia R package85.