Ecologists collectively predict that climate change will enhance phytoplankton biomass in northern lakes. Yet there are unique variations in the structures and regulating functions of lakes to make this prediction challengeable and, perhaps, inaccurate. We used archived Landsat TM/ETM+ satellite products to estimate epilimnetic chlorophyll-a concentration (Chl-a) as a proxy for phytoplankton biomass in 281 northern temperate lakes over 28 years. We explored the influence of climate (air temperature, precipitation) and landscape proxies for nutrient sources (proportion of wetlands in a contributing catchment, size of the littoral zone, potential for wind-driven sediment resuspension as estimated by the dynamic ratio) or nutrient sinks (lake volume) in a random forest model to explain heterogeneity in peak Chl-a. Lakes with higher Chl-a (median Chl-a = 2.4 μg L^-1, n = 40) had smaller volumes (< 44 × 10⁴ m³) and were more sensitive to increases in temperature. In contrast, lakes with lower Chl-a (median Chl-a = 0.6 μg L^-1, n = 241) had larger volumes (≥ 44 × 10⁴ m³), contributing catchments with smaller proportions of wetlands (< 4.5% of catchment area, n = 70), smaller littoral zones (< 16.4 ha, n = 137), minimal wind-driven sediment resuspension (as defined by the dynamic ratio; < 0.45, n = 232), and were more sensitive to increases in precipitation. Lakes with larger volumes were generally less responsive to climate factors; however, large volume lakes with a significant proportion of wetlands and larger littoral zones behaved similarly to lakes with smaller volumes. Our finding that lakes with different landscape properties respond differently to climate factors may help predict the susceptibility of lakes to eutrophication under changing climate conditions.

Air temperature and precipitation data were extracted from 60-arcsecond historical monthly air temperature and precipitation grids for each of 281 study lakes. Landscape factors (e.g., lake volume, depth, lake catchment size) were extracted using Digital Terrain Models (DEMs) for the province of Ontario (Canada) and georeferenced and digitized countour maps. The calculation of all factors was performed in ArcGIS 10.8 and Excel. Statistical analysis was performed in R 3.6.0.