Habitat-Mediated Responses of Zooplankton to Decreasing Light in Two Temperate Lakes Undergoing Long-Term Browning

Increases in dissolved organic matter and the consequent “browning” of some lakes in recent decades are reducing water transparency to both ultraviolet and photosynthetically active radiation with important, but poorly understood ecosystem-level consequences for zooplankton grazers. The prevailing resource-based unimodal hypothesis posits that nutrients in dissolved organic matter stimulate primary production in clear-water lakes, while shading by dissolved organic matter inhibits primary production in browner lakes, with zooplankton responses following the patterns of their food resources. Support for this hypothesis derives primarily from short-term experiments, space-for-time analyses, and modeling studies. Here we use three decades of long-term monitoring data from two temperate lakes to assess zooplankton responses to changes in not only resources (chlorophyll) as drivers of change, but also light-related habitat variables (ultraviolet and photosynthetically active radiation, surface and deep-water temperatures, deep-water dissolved oxygen, and pH). The study lakes include one clear-water lake and one browner lake, both of which have experienced long-term browning. Given that zooplankton depth distribution and body size can vary with water transparency, color, and temperature, we test for responses in not only overall zooplankton abundance, but also in vertical distribution and body size. We also examine the ability of the relationship between short-term interannual variation in the driver vs. response variables to predict long-term zooplankton trends. The primary responses of zooplankton were strong changes in abundance that varied with taxon and life history stage in response to habitat variables rather than food resources. Only two groups showed vertical distribution responses, and they trended toward deeper distributions with browning. There was no significant change in body size ratio. The directionality of the response of zooplankton abundance to interannual variability in the driver variables was consistent with those observed in the long-term trends for 33 of 80 comparisons (41%), though only three of those (4%) had statistically significant short-term interannual variability relationships. We conclude that habitat-related changes associated with browning in lakes have important consequences for zooplankton community structure with stronger effects in clear-water lakes than in browner lakes, and that even at the whole-lake scale short-term data are generally not adequate to predict long-term responses.

Increases in dissolved organic matter and the consequent "browning" of some lakes in recent decades are reducing water transparency to both ultraviolet and photosynthetically active radiation with important, but poorly understood ecosystem-level consequences for zooplankton grazers. The prevailing resource-based unimodal hypothesis posits that nutrients in dissolved organic matter stimulate primary production in clear-water lakes, while shading by dissolved organic matter inhibits primary production in browner lakes, with zooplankton responses following the patterns of their food resources. Support for this hypothesis derives primarily from short-term experiments, space-for-time analyses, and modeling studies. Here we use three decades of long-term monitoring data from two temperate lakes to assess zooplankton responses to changes in not only resources (chlorophyll) as drivers of change, but also light-related habitat variables (ultraviolet and photosynthetically active radiation, surface and deep-water temperatures, deep-water dissolved oxygen, and pH). The study lakes include one clear-water lake and one browner lake, both of which have experienced long-term browning. Given that zooplankton depth distribution and body size can vary with water transparency, color, and temperature, we test for responses in not only overall zooplankton abundance, but also in vertical distribution and body size. We also examine the ability of the relationship between short-term interannual variation in the driver vs. response variables to predict long-term zooplankton trends. The primary responses of zooplankton were strong changes in abundance that varied with taxon and life history stage in response to habitat variables rather than food resources. Only two groups showed vertical distribution responses, and they trended toward deeper distributions with browning. There was no significant change in body size ratio. The directionality of the response of zooplankton abundance to interannual variability in the driver variables was consistent with those observed in the long-term trends for 33 of 80 comparisons (41%), though only three of those (4%) had statistically significant short-term interannual variability relationships. We conclude that habitat-related changes

INTRODUCTION
Water clarity is changing in lakes worldwide due to anthropogenic and natural environmental factors ranging from lake ontogeny following glacial recession (Engstrom et al., 2000;Williamson et al., 2001;Sommaruga, 2015) to increases in eutrophication and browning that are related to extreme precipitation events, changes in land use, and other anthropogenic activities (Kritzberg and Ekström, 2012;Williamson et al., 2016Williamson et al., , 2017Kritzberg, 2017;Rose et al., 2017;Leech et al., 2018). One of the most widespread changes in water clarity is the browning of lakes due to increases in the quantity and quality (color) of terrestrially-derived dissolved organic matter (DOM, often measured as concentration of dissolved organic carbon, DOC) (Monteith et al., 2007;SanClements et al., 2012;Kritzberg, 2017;Leach et al., 2019). This browning is expected to continue into the future with climate change-related increases in precipitation and extreme storm events (Larsen et al., 2011;de Wit et al., 2016;Weyhenmeyer et al., 2016).
Two largely light-driven hypotheses, resource-, vs. habitatbased, have been used to explain changes in zooplankton consumers in response to browning (Figure 1). The resourcebased, unimodal hypothesis argues that during browning, basal resources are the primary driver of changes in phytoplankton Seekell et al., 2015a,b), zooplankton (Jones et al., 2012;Kelly et al., 2014Kelly et al., , 2016, and fish (Finstad et al., 2014). According to this hypothesis, in clear-water lakes with low DOM (generally < ∼5 mg C L −1 ), increases in DOM introduce nutrients that stimulate increases in primary production and consequently zooplankton production, while in lakes above this threshold the effects of DOM shading of PAR reduce primary production and consequently zooplankton production. Support for this hypothesis derives from short-term experiments, spacefor-time analyses, and modeling .
While nutrients and light control primary production in the resource-based hypothesis, in the habitat-based hypothesis, light alters key vertical habitat characteristics such as temperature, oxygen, and the potential for UVR damage. In addition to the direct reductions of UVR and PAR, declines in water transparency to PAR result in warmer surface waters, cooler temperatures and depleted dissolved oxygen in deep waters, and pH also increases, all of which are typical of lakes undergoing browning in northeastern North America and northern Europe (Monteith et al., 2007;Williamson et al., 2015;Strock et al., 2017). One study found that benthic consumers are more influenced by habitat changes such as oxygen depletion than by changes in resources in lakes spanning a range of DOC concentrations (Craig et al., 2015). Similarly, an analysis of ∼20 years of data on a suite of 28 Adirondack lakes that have been undergoing browning reported that changes in pelagic primary producers were more related to changes in light than to nutrients, and the responses of phytoplankton and zooplankton to browning differed from each other, reflecting a trophic decoupling rather than a resource-based response of the zooplankton directly to changes in phytoplankton (Leach et al., 2019). These studies extend earlier work that has shown the importance of light limitation of primary and secondary production with increasing DOM Kelly et al., 2014;Thrane et al., 2014).
While it is widely recognized that browning decreases both photosynthetically active (PAR) and ultraviolet (UVR) radiation in lakes (Carpenter et al., 1998;Finstad et al., 2014;Kelly et al., 2014;Seekell et al., 2015a), the effects of decreasing UVR are rarely considered due to the lack of available instrumentation and data. Yet decreases in UVR exposure can be very pronounced during browning due to the selective absorption of shorter wavelength UVR by terrestrially-derived DOM. In seepage lakes with minimal riverine inputs, DOM is the primary regulator of water transparency to UVR, and also to PAR (Morris et al., 1995). In addition, the non-linear relationship between DOC concentration and the compensation depth (1% of subsurface PAR, the depth above which there is net photosynthesis, and below which there is net respiration) or equivalent 1% depth to which UVR penetrates (Williamson et al., 1996) causes clear, low DOC lakes to experience much stronger changes in their underwater light environment during browning than lakes with initially higher DOC concentrations (Snucins and Gunn, 2000;Read and Rose, 2013;Williamson et al., 2015). Thus in clearwater lakes even small increases in DOM strongly reduce the underwater UVR exposure levels (Williamson et al., 1996), with consequences that have the potential to alter zooplankton nutrition (Nova et al., 2019), predation (Williamson et al., 1999;Boeing et al., 2004;Leech et al., 2009;Lindholm et al., 2016), parasite loads (Overholt et al., 2012;Williamson et al., 2017), behavior Wolf and Heuschele, 2018), and vertical distribution in the water column (Leech and Williamson, 2001;Cooke et al., 2008;Fischer et al., 2015;Leach et al., 2015;Urmy et al., 2016). Here we have a unique, long-term database that includes not only well-resolved changes in PAR during browning, but also changes in UVR. This enables us to test whether substantial decreases in both PAR and UVR are, in fact, related to changes in zooplankton distribution and abundance.
Our core objective here is to investigate the effects of changes in the underwater light environment, including UVR as well as PAR, on zooplankton consumers in lakes undergoing browning. We use long-term data from Lake Giles and Lake Lacawac, two browning temperate lakes in northeastern Pennsylvania, USA, FIGURE 1 | Conceptual model of how DOM alters the underwater light environment to influence not only resources relevant to the resource-based hypothesis (green), but also many aspects of the pelagic habitat for zooplankton (black). Direct effects of PAR and UVR (e.g., vertical migration, DNA damage) are shown in red.
to test the importance of light-related drivers on zooplankton, and assess the relative importance of resource-based vs. habitatbased responses in the pelagic zone of lakes (Figure 1). Giles and Lacawac are both seepage lakes that have experienced longterm changes in their underwater light environments, with Giles responding more strongly to browning than Lacawac due to its initially very low DOC concentrations . Here, we extend this prior study by (1) adding 5 years to the longterm data set, (2) analyzing the zooplankton data with higher taxonomic and life stage resolution, as well as with additional metrics including body size and vertical distribution since these variables have been shown to vary with water transparency, water temperature, and water color, and (3) relating changes in zooplankton to both the interannual variability (IAV) and the long-term trends (LTT) in light-related driver variables, following the approach of Leach et al. (2019). Importantly we examine whether using IAV can provide critical early insights into population-and ecosystem-level responses of zooplankton to both abiotic and biotic forcing over the long term. While short-term studies provide valuable insights into changes in lake structure and function in response to disturbance or extreme events, their ability to predict long-term responses of consumers is less certain.

Sampling Methods
Lakes Giles and Lacawac are two small lakes in northeastern Pennsylvania, USA with surface areas of 48.0 and 21.4 hectares, maximum depths of 24 and 13 m, and mean depths of 10.1 and 5.2 m, respectively. Water samples were collected for pH, DOC, chlorophyll, and nutrient analyses from the epilimnion of both lakes using a Van Dorn sampler. Data from June, July and August were averaged to provide a single value for each year from 1988 to 2019. For measurements of DOC concentration, samples were filtered through a pre-ashed 0.7 µm glass fiber filter and analyzed with a TOC analyzer using a high temperature oxidation method. Water samples were analyzed for pH using an Orion pH meter, calibrated at 4.00 and 7.00 using commercial pH buffers. For measurements of total phosphorus, unfiltered samples were acidified using H 2 SO 4 , digested with a potassium persulfate solution, and analyzed using the molybdenum blue method. For measurements of soluble reactive phosphorus (SRP), samples were filtered through an A/E glass fiber filter (1 µm), acidified with H 2 SO 4 and analyzed using the molybdenum blue method. Chlorophyll samples were filtered onto 1.0 µm A/E filters and immediately frozen. Extraction methods for chlorophyll varied among years, but followed standard methodologies as detailed in Williamson et al. (2015). Using a YSI temperature-oxygen meter (1988)(1989)(1990)(1991)(1992) or a submersible temperature and UVR-PAR radiometer (1993-2019), epilimnetic temperature was measured at a depth of 2 m in both lakes, while hypolimnetic temperature was measured at depths of 18 m (Giles) and 10 m (Lacawac). The YSI temperature-oxygen meter was also used to measure dissolved oxygen (all years) at 1 m or greater resolution. The deep-water percent saturation of dissolved oxygen (20 m in Giles, 11 m in Lacawac) was calculated from dissolved oxygen concentration and water temperature using the "rMR" R package (Moulton, 2018). Light profiles were collected using submersible radiometers to measure attenuation of UVR (µW cm −2 nm −1 at 320 nm, with full width at half maximum bandwidth = 8 nm) and PAR (µmol m −2 s −1 400-700 nm). The radiometers used included Biospherical Instruments (BSI, San Diego, California) PUV radiometers (PAR and UVR;1993-2003 and BSI BIC radiometers with a deck cell (PAR and UVR;2003. The depth to which 1% of subsurface irradiance penetrated was derived from the vertical profile data for PAR and UVR using the diffuse attenuation coefficient, K d , calculated as the slope of the log-linear relationship between irradiance and depth (Williamson et al., 1996).
Zooplankton were collected according to the methods detailed in Williamson et al. (2015), and daytime densities from June, July, and August were averaged to provide a single annual value. Briefly, net tows with a closing bongo net (side-byside Wisconsin-style 48 µm and standard 202 µm mesh nets) encompassed the full water column at the deepest point in each lake. Most often zooplankton samples were collected separately from the epilimnion and the meta-hypolimnion using the net's closing feature. These two strata were distinguished based on thermal profiles where the bottom of the epilimnion was defined as where temperature change was >1 • C per m depth. For long-term trend analysis, samples from all depths were combined using depth-weighted densities to provide a single water column value. Daphnia from either a 48 µm or a 202 µm net were counted in a Bogorov chamber. Comparison of Daphnia abundance using samples from the 48 and 202 µm counts showed no difference in estimates of abundance. Calanoid and cyclopoid copepods as well as rotifers were collected with a 48 µm mesh net and counted in a Sedgewick-Rafter chamber. Copepod adults and copepodids (juveniles) were separated into cyclopoids or calanoids. Cyclopoid and calanoid nauplii were counted together as nauplii and not differentiated. Total calanoids and total cyclopoids included adults and copepodids of the respective taxon, but not nauplii. The proportion of zooplankton in the epilimnion was estimated by dividing the proportion of a given zooplankton taxon collected in the epilimnion by the proportion of the depth of the water column that comprised the epilimnion. A value near 0.5 indicates an approximately even distribution of the zooplankton taxon between the epilimnion vs. deeper depths (metalimnion plus hypolimnion), while a larger value (e.g., nearer to 1) indicates a more epilimnetically-distributed taxon, and a smaller value (e.g., nearer to 0) indicates a more deeply distributed taxon. Zooplankton body size ratios (large: small) were calculated as the sum of all Daphnia and adult cyclopoid and calanoid copepods divided by the sum of all rotifers and nauplii.

Statistical Analysis
Following the data collection and preparation described above, the eight primary driver variables incorporated into the analyses had a median of 27 years of data out of a maximum of 32 (1988 through 2019), with a minimum of 19 years of data to a maximum of 32 years of data (Supplementary Table 1). The zooplankton response variables tended to have fewer years of data, with a median of 17 years of data, a minimum of 12 and a maximum of 20 (Supplementary Table 1). A nonparametric Mann-Kendall test was used to determine if there were statistically significant trends over time for each driver and zooplankton response variable for each lake using the full set of annually-averaged data with α = 0.05, and Sen's slope was calculated to determine the direction and rate of change over time, using the "trend" R package (Pohlert, 2018). Mann-Kendall non-parametric trend tests and Sen's slope estimator have relaxed assumptions of normality and linearity compared to parametric tests, allow missing data points, and are robust to outliers and data gaps (Gilbert, 1987). Interannual variability correlations were conducted following the methods of Leach et al. (2019). Briefly, the first derivative of each variable was taken over time. Data gaps >1 year were not included in the subsequent correlations. Non-parametric Spearman rank correlations were conducted on the first derivative data for each pair of driver and zooplankton response variables independently for each lake. The resulting correlation coefficient (ρ) for IAV was compared to the pair of LTT to determine if there were consistent patterns at the interannual vs. long-term time scales. Consistent patterns included when long-term Sen's slope trends for the pair of abiotic and zooplankton LTT responses had the same sign (+/+ or −/−) and a positive ρ for IAV, or when long-term Sen's slope trends for the pair had opposing signs (+/−) and a negative ρ for IAV. If any given pair of driver and response variables had a significant LTT with a consistent directional pairing of LTT slopes and IAV ρ, it was considered to have some predictive power to estimate LTT from IAV. If these conditions were met and the IAV ρ was statistically significant (α = 0.05), it was considered to have stronger predictive power (a higher level of confidence that IAV can predict LTT). All figures were created using the "ggplot2" R FIGURE 2 | Long-term summer average (June, July, August) trends of PAR (A) and 320 nm UVR radiation (B), DOC (C), and chlorophyll a (D) in Lake Giles (blue, open circles with solid line), and Lake Lacawac (brown, filled circles with dashed line). There was an abrupt change in PAR in Lake Giles following a strong variation in 2003-2004 that provided contrasting underwater PAR environments in the periods before 2003 and after 2004. Trends of change in UVR were much more gradual, and not as strong in Lacawac. Letters indicate statistically significant long-term trends in Giles ("G") and in Lacawac ("L"; Mann-Kendall non-parametric trend test, p < 0.05). Black trend lines are LOESS smoothed lines.
package (Wickham, 2016), and all analyses were completed in R version 3.6.2 (R Core Team, 2019).

RESULTS
Adding 5 years to our previous data set shows that over the past three decades the LTT in compensation depth (1% PAR depth) in Lake Giles consisted of a decrease from 20 m or more to <15 m on average (Figure 2A). The 1% 320 nm UVR depth shallowed from close to 10 m to around 2.5 m ( Figure 2B). Water clarity was high from 1993 to 2002, highly variable from 2003 to 2004, and substantially lower from 2005 to 2019 (Figure 2A). Lacawac had decreases in both PAR and UVR penetration that were much less pronounced than those in Giles, with only the trends in UVR being statistically significant (Figure 2, Supplementary Table 1).
There have been strong and significant increases in the LTT in DOC concentration in Lake Giles (Figure 2C), where DOC has more than doubled from ∼0.9 to 2.5 mg C L −1 over three decades. Increases in DOC were much less pronounced in Lacawac and not statistically significant ( Figure 2C). While DOC has increased, chlorophyll has been highly variable with no consistent change over time in either lake ( Figure 2D). Even FIGURE 3 | Vertical profiles of PAR (A) and 320 nm UVR (B) in Lake Giles relative to the depth of the mixed layer (dashed lines) in 1993 and 2019. Note the strong shallowing of the compensation depth (1% of subsurface) for PAR, which influences the balance of photosynthesis and respiration, and thus the depth of light limitation of photosynthesis as well as oxygen depletion in the deeper waters. There is also a strong reduction in UVR and PAR exposure levels in the mixed layer.
at deeper depths, where we've observed a deep chlorophyll maximum in the past, the only significant trend was a decrease in chlorophyll in the hypolimnion of Lake Giles with browning (Supplementary Table 1). Similar tests on the 10-12 years of nutrient data indicated no significant change in total or soluble reactive phosphorus (Supplementary Table 1). The DOC: TP ratio did show a significant increase in Lake Giles, but not in Lacawac, consistent with a significant increase in DOC in Giles but not Lacawac, and no trend in TP in either lake (Supplementary Table 1).
In contrast to the lack of change in chlorophyll and nutrients, there were critical light-driven changes in habitat variables during browning, especially in Lake Giles. Lower water transparency reduced both sunlight penetration and surface mixed layer depth. For example, early in the study period (1993) in Lake Giles the compensation depth (1% PAR) extended to the bottom of the lake, supporting net photosynthesis throughout the water column. In contrast, after nearly three decades of browning (2019) the 1% PAR depth reached only about half as deep as in 1993 (Figure 3A). Similarly, in 1993, potentially damaging 320 nm UVR was 10% of the subsurface irradiance or greater throughout the surface mixed layer (Figure 3B), but by 2019, the mixed layer itself was only about half as deep, while the 1% UVR depth decreased by more than four-fold ( Figure 3B). Corresponding to the strong decreases in PAR in Giles between 2002 and 2005 were particularly strong decreases in deepwater temperature and dissolved oxygen, as well as increases in surface water temperature and pH (Figures 2A, 4). These primarily light-driven responses led to important ecosystem-level consequences including a significant decrease in chlorophyll in the hypolimnion (Supplementary Table 1), the onset of oxygen depletion in the hypolimnion ( Figure 4C) and a refuge from damaging UVR exposure in the deeper portions of the warmer surface waters (Figure 3B). Lacawac showed a variable but significant decline in hypolimnetic temperature but no trends in surface temperature, dissolved oxygen, or pH (Figure 4).
Zooplankton LTT showed a diversity of changes in abundance, some changes in their daytime vertical distribution, but no changes in the community body size ratio (Supplementary Table 1). Similar to the changes in habitat variables, zooplankton responses were stronger in Giles than in Lacawac. Zooplankton were generally most abundant in the deeper waters, or distributed approximately evenly throughout the water column in both lakes ( Table 1). Only the rotifer Keratella was distinctly epilimnetic (Table 1), and only two groups showed statistically significant changes in their vertical distribution (Figure 5). Daphnia in Lacawac showed a statistically significant but highly variable decrease in their proportion in the epilimnion (p = 0.034, Figure 5A, Supplementary Table 1), and nauplii in Giles showed a stronger and more consistent decrease in their proportion in the epilimnion (p = 0.004, Figure 5B, Supplementary Table 1). Copepods, cladocerans, and rotifers all exhibited significant decreasing (calanoids, Daphnia, Gastropus, Polyarthra) or increasing (cyclopoids, Kellicottia) LTT in abundance in Lake Giles (Figure 6). In Lacawac, abundance changes were only observed in Gastropus, which also declined ( Figure 6G), and in cyclopoid copepod adults, which increased significantly ( Figure 6D). The rotifer Keratella showed no significant LTT in abundance in either lake (Supplementary Table 1).
Of the 80 pairings that compared zooplankton responses to each potential driver variable on both an interannual vs. longterm basis, 33 (41%) showed consistent IAV and LTT responses (e.g., both positive, both negative, or both opposite). Only three (4%) of the total showed consistent IAV and LTT responses that were also statistically significant, all in Giles (Figure 7). Two of these three strong relationships involved light as a driver variable (one UVR and one PAR), while one involved pH (Figure 7). Chlorophyll a, an estimate of resource availability, was the only variable with no explanatory power, and no consistent relationships between IAV and LTT, due primarily to the lack of any significant LTT.  (Pilla et al., 2018). The increases in pH (D), statistically significant only in Lake Giles, are most likely due to decreases in acid deposition in the region during the study period . Letters indicate statistically significant long-term trends in Giles ("G") and in Lacawac ("L"; Mann-Kendall non-parametric trend test, p < 0.05). Black trend lines are LOESS smoothed lines.

DISCUSSION
Browning has induced strong changes in the underwater light environment in lakes that have in turn altered pelagic habitats with important consequences for zooplankton consumers. The light-related changes in habitat and responses by zooplankton observed here were much stronger than any changes in basal resources (i.e., chlorophyll), and also much more prevalent in clear-water Lake Giles than in browner Lake Lacawac. The resource-based hypothesis argues that nutrients in terrestrially-derived DOM are responsible for increases in primary production in clear-water lakes, while in browner lakes shading plays a more important role, and leads to a decrease in primary production. This leads to a unimodal relationship where in clear-water lakes such as our Lake Giles, and possibly more intermediate productivity Lake Lacawac, one would predict that as DOC concentration increases, nutrients and thus chlorophyll would also increase. Yet phosphorus did not change in either of the study lakes, and the only significant trend in chlorophyll was a decrease in the hypolimnion of Lake Giles with browning, which is the opposite of what would be expected according to the resource-based hypothesis. Chlorophyll is only a proxy for food resources, and the chlorophyll per unit biomass can increase at lower light levels. The fact that chlorophyll exhibited no increases even with browning-related reductions in light is thus also in contrast to the predictions of the resource-based hypothesis. Another possible, often mentioned but rarely tested, resourcebased control of primary production at low levels of DOM is photoinhibition by UVR (Carpenter et al., 1998;Finstad et al., 2014;Kelly et al., 2014;Seekell et al., 2015a). Photoinhibition, which has been demonstrated in situ in Giles (Moeller, 1994), should be reduced with increases in DOM, leading to increases in primary production. However, the lack of any significant increases in chlorophyll or nutrients suggest that UVR effects on basal resources are also not the most likely explanation for the observed changes in zooplankton in these two lakes. Changes in vertical habitat gradients were more important.
We observed LTT in both abiotic driver variables and zooplankton responses. In only 4% of the cases, however, were LTT trends highly predictable from IAV, making it challenging to use short-term data to predict long-term responses. The ability to assess LTT based on short-term experiments or IAV in whole lakes could provide early warning signals that enable implementation of management options to avoid the most severe changes in ecosystems well before they occur (Pace et al., 2015). Analyses that use a single or subset of only a few extreme years (e.g., heat waves or precipitation extremes) provide a convenient natural manipulation that promotes insights into how environmental change can alter lake ecosystem structure and function. For example, a single extremely warm year during a heat wave in Europe showed much stronger thermal stratification of the water column and deep-water oxygen depletion than normal (Jankowski et al., 2006). Comparison of 86 small Canadian lakes in one cool and two very warm years showed changes in thermal structure that were strongest in the clearer lakes (Snucins and Gunn, 2000). A study of 84 lakes in northeastern North America showed that an extreme wet year can lead to higher DOC concentrations and short-term browning, while an extreme dry year can lead to higher sulfate concentrations and acidification (Strock et al., 2016). Wholelake manipulations can provide a similar short-term snapshot of how lake ecosystems respond to short-term changes (Christensen et al., 1996), though longer-term changes may not be detectable for many decades (Carpenter and Pace, 2018). The presence of so few consistent relationships between zooplankton responses and potential drivers at these two different time scales of IAV and LTT for the majority of our analyses (59%), however, suggests caution is needed before extrapolating LTT from short-term experiments or data, even at the whole-lake scale. We did observe both direct and indirect effects of a changing light environment on pelagic habitats during browning, and they are likely related to regional increases in precipitation in addition to decreases in acid deposition. While there has been no significant LTT in air temperature or solar radiation in this region during the study period, there has been an increase in average precipitation with a shift to a wetter regime between 2001 and 2004 (Pilla et al., 2018), a period of particularly rapid change in light attenuation. In addition to the increases in DOC concentration documented here, there has been a precipitation-related increase in the color (DOC-specific UVR absorbance) of the DOM (Williamson et al., 2014). Thus, the changes in the underwater light environment observed here are a result of not only the observed increases in DOC concentration, but also increases in color.
Our unique long-term UVR database enabled us to examine the role of UVR-related changes and responses during browning. UVR has the potential to have direct negative effects on zooplankton. A recent meta-analysis on the negative effects of UVR on a wide variety of freshwater organisms ranging from phytoplankton and zooplankton to amphibians and fish, showed the greatest negative effects on zooplankton (Peng et al., 2016). High exposure levels of >10% of surface 320 nm UVR, such as those observed historically in Lake Giles (Figures 2B, 3), have been related to the elimination of all zooplankton species in shallow montane lakes in Argentina, with the exception of one very highly pigmented copepod (Marinone et al., 2006). This high UVR exposure level has been experimentally demonstrated to FIGURE 7 | Relationships between IAV and LTT of driver variables vs. zooplankton response variables. Rows are potential driver variables related to the resource-based or habitat-based hypotheses. Columns are zooplankton abundance response variables that had a significant long-term trend in either Giles or Lacawac. Small circles indicate when both the respective driver variable and zooplankton response variable LTT were statistically significant (Mann-Kendall non-parametric trend test, p < 0.05), and the direction of the interannual Spearman correlation was consistent with that of the LTT for the same lake (open white circles = Giles, filled black circles = Lacawac). Large points follow the same pattern, but indicate when the interannual Spearman correlation was consistent with that of the LTT, and also statistically significant (Spearman rank correlation, p < 0.05), suggesting a higher degree of confidence in predicting LTT from IAV data. Signs indicate the positive (+) or negative (-) nature of the relationship between the driver and response variables.
lead to high mortality rates of a range of not only zooplankton (Williamson et al., 1994(Williamson et al., , 1999, but also early life history stages of fish (Williamson et al., 1999;Huff et al., 2004;Olson et al., 2006Olson et al., , 2008. However, the decreases in Daphnia, Polyarthra, and Gastropus during a period of decreasing UVR exposure in the surface waters are not consistent with a direct UVR-damage hypothesis (Leech and Williamson, 2000;Persaud and Williamson, 2005). Similarly, zooplankton may also have historically been squeezed between UVR in the surface waters, and deeper-dwelling tactile invertebrate predators (Boeing et al., 2004), but we would expect a release from this squeeze with declining UVR, which contrasts with the observed declines in abundance of zooplankton grazers. Zooplankton have a wide variety of mechanisms to reduce UVR damage ranging from behavioral avoidance to the production of photoprotective compounds and photoenzymatic repair, though these protections come at some cost (Rautio and Tartarotti, 2010). While few zooplankton inhabited the surface waters of either of our lakes (Table 1), no increases in abundance in the epilimnion were observed with the decreases in UVR over time as might be expected if direct UVR damage was an important driver. Furthermore, the rotifer Keratella was the most epilimneticallydistributed of all of the species, and is also known to be one of the most UVR-tolerant zooplankton species in Giles (Williamson et al., 1994). Keratella showed no significant LTT in abundance and no change in vertical distribution with browning in either lake. In spite of the strong declines in surface UVR, it was the deep-water taxa such as cyclopoids and Kellicottia that showed the greatest increases over time, though they were not exposed directly to high levels of UVR at the depths where they were found. In contrast, calanoid copepods from Lake Giles are attracted to  or benefit from (Cooke and Williamson, 2006) exposure to UVR. This is consistent with declines in calanoid copepods observed with concurrent declines in UVR exposure in Giles. In both Giles and Lacawac multiple single or interacting habitat-related factors are likely more important to the response of zooplankton than is direct exposure to UVR damage.
There is a strong interplay between reductions in water transparency to PAR, increases in surface water temperatures, and decreases in deep-water temperatures and oxygen. For example, a multi-lake modeling study, which included Giles, showed that in lakes where both transparency and mixed layer depth were decreasing, the overall exposure to UVR and PAR in the mixed layer decreased, whereas in lakes where only the mixed layer depth decreased, average exposure increased (Neale and Smyth, 2019). Coupling these changes in lightrelated habitat variables to what is known about the ecology of zooplankton has the potential to add insights to the changes in zooplankton abundance and depth distribution presented here. Physiological tolerances to temperature and low oxygen involving cyclopoid copepods provide a good example of the importance of these complex interactions. Cyclops scutifer is the dominant cyclopoid in Giles, and it has shown some of the strongest increases with browning of any of the zooplankton taxa (Figures 6C,E). This may be related to several factors. First, C. scutifer is a cold stenotherm that inhabits deep waters during warmer seasons (Elgmork, 1967;Johnson et al., 2007). The rotifer Kellicottia can increase under colder conditions as well (Persaud and Williamson, 2005). One potential threat to cold stenotherms is being squeezed between two suboptimal habitats with warming surface waters, and oxygen-depleted deeper waters (Doubek et al., 2018), both of which have been observed in Lake Giles (Figures 4A,C). However, the deep waters of Giles have been cooling (Figure 4B), and oxygen depletion has not yet reached a state of anoxia as it has in Lacawac ( Figure 4C). Therefore, the increase in optimal deep, cold-water habitat may actually be contributing to the increase in both Kellicottia and C. scutifer. The deeper depth distribution of nauplii with browning in Giles ( Figure 5B) is most likely related to a shift from calanoid nauplii to cyclopoid nauplii as hypolimnetic cyclopoids increased in abundance and the calanoids declined ( Figure 6, Table 1). Some cyclopoids such as Mesocyclops edax can tolerate anoxic conditions for at least short periods of time, and even use anoxic waters as a refuge from fish predation (Williamson and Magnien, 1982). M. edax is the species that is most abundant in the deeper hypoxic to anoxic waters of Lacawac ( Figure 4C, Table 1). If hypoxia continues to intensify in Giles in the future, we may see a shift to this more anoxiatolerant species.
While recovery from acidification and an increase in pH are characteristic of browning in lakes (Evans et al., 2006;Monteith et al., 2007;Williamson et al., 2015;Strock et al., 2016), shortterm experiments that manipulate DOM from the microcosm to the whole-lake level as well as comparative space-for-time studies of browning have rarely considered changes in pH. Yet there is good evidence that zooplankton species composition varies with pH. A comparative space-for-time study of 20 lakes in northeastern North America showed a loss of biodiversity with a decrease in pH (Confer et al., 1983). While Leptodiaptomus minutus was found in acidified lakes, other zooplankton species such as Daphnia and cyclopoid copepods were rarely if ever found at a pH below 5. A comparative study of 32 lakes in Glacier Bay, Alaska, with a pH range from 5 to 8 found Cyclops scutifer and rotifers to be most abundant in lakes with a higher pH (Olson et al., 1995). During an experimental whole-lake basin acidification experiment, Gastropus stylifer was found to be most abundant at the most acidic pH of 4.7, while Kellicottia longispina was found at greater abundances at the two higher pH values of 5.2 and 5.6 (Frost et al., 1998). These patterns are consistent with the observed increases over time in cyclopoid copepods in Giles and Lacawac (Figures 6C,E) as well as Kellicottia in Giles (Figures 6H, 7) where pH has also been increasing, though not significantly in Lacawac ( Figure 4D). Increases in pH likely also contributed to the increases in fish that have been observed in these study lakes , as fish tend to increase above a pH of around 5 in lakes recovering from acid deposition (Baldigo et al., 2016).
A third hypothesis that may contribute to changes in zooplankton communities during browning that has received less attention, but for which good evidence is accumulating, is predation. The three most important zooplanktivorous predators in our study lakes are cyclopoid copepods, larvae of the phantom midge Chaoborus, and young-of-year fish. Only our data on cyclopoid predators was sufficient to incorporate into the full analysis here, though we also have very limited data on young-of-year fish. Below we discuss how predation by both tactile invertebrates as well as visual predators has the potential to be an important indirect effect of browning on zooplankton.
Prior experimental work suggests that predation regimes will shift from being dominated by visual predators to domination by tactile invertebrate predators with increases in DOM (Wissel et al., 2003a). While the increase in abundance of cyclopoids in Giles is consistent with this prior experimental work, we also observed an increase in visual predators in our lakes . As discussed above, these increases in cyclopoids and young-of-year fish may be due at least in part to the mitigating effects of recovery from acid deposition and a higher pH. This highlights the importance of browningrelated factors other than resources or light alone in mediating the observed changes in zooplankton, including indirect effects such as predation.
The effects of predatory cyclopoids on zooplankton are most likely to be on smaller zooplankton species, and the rotifers in particular. Rotifers are an important food source for many cyclopoid copepods (Williamson, 1983(Williamson, , 1984. Cyclopoid copepod predators in lakes with warming surface waters (as observed in our lakes) can also benefit from increased availability of rotifers as food . The changes that we observed in rotifer species are consistent with increases in cyclopoid predation. For example, the spined rotifers Keratella and Kellicottia are well-protected against cyclopoid predation (Figure 8), and either did not change in abundance, or increased, respectively ( Figure 6H). In contrast, the soft-bodied rotifers Gastropus and Polyarthra, which are poorly defended against cyclopoid predators (Gilbert and Williamson, 1978; Figure 8), declined (Figures 6G,I).
Larvae of the phantom midge Chaoborus are some of the most important invertebrate predators of Daphnia in lakes (Riessen et al., 2012). Chaoborus are, however, highly heterogeneous in their horizontal distribution and migrate into the sediments during the day. Thus, adequate sampling requires multi-station sampling at night (Wissel et al., 2003b), something which was beyond our capacity in this long-term study. There is evidence, however, that Chaoborus abundance FIGURE 8 | Photographs showing a cyclopoid copepod (Mesocyclops edax) (A) and differential defenses of rotifers against cyclopoid predation including the spined rotifers Kellicottia (B) and Keratella (C) that either increased or did not change in abundance, respectively, and the soft-bodied rotifers Gastropus (D) and Polyarthra (E) that showed significant declines. Size of rotifers magnified relative to the copepod to show detail.
is positively correlated with DOM (Wissel et al., 2003a,b), and thus an increase in Chaoborus with browning would be expected. A recent paper has provided evidence that the refuge from UVR due to increased DOM enables Chaoborus in shallow arctic ponds to decimate their zooplankton prey (Lindholm et al., 2016). Thus, increases in Chaoborus are one possible explanation for the observed decline in Daphnia in Giles.
Changes in the underwater light environment can also influence visual predation by fish on larger zooplankton such as Daphnia. While increasing DOM can induce light limitation that reduces predation by bluegill (Lepomis macrochirus) and largemouth bass (Micropterus salmoides) (Weidel et al., 2017), there are more factors influencing predators than just visible light limitation during long-term lake browning. For example, bluegill and largemouth bass are two of the dominant fish predators in our lakes, yet sampling of young-of-year fish early and later during the process of browning has shown dramatic increases in fish abundance in both of our lakes, with an increase in catch-per-unit effort of young-of-year fish from 0 and 1.6 in 1991 in Giles and Lacawac respectively, to 5.0 and 5.5 in 2014 . Early life history stages of planktivorous fish have historically been highly susceptible to UVR damage in Lake Giles (Williamson et al., 1997(Williamson et al., , 1999Huff et al., 2004), and have increased in the past three decades as UVR exposure levels have declined . This increase in fish may also contribute to the observed increase in cyclopoids because cyclopoids are known to be less susceptible to visual predators than are cladocerans, and thus are more abundant in lakes with fish (Johnson et al., 2010). The hypolimnetic nature of the cyclopoids in our lakes (Table 1) would additionally make them less susceptible to visual predators like young-of-year fish that primarily inhabit the warmer surface waters. These patterns support a stronger habitat-driven rather than resource-driven effect of browning when the predation changes and taxonspecific and life history stage-specific changes in zooplankton are considered.
Although the patterns of change in zooplankton observed during browning in this study are not consistent with the resource-based unimodal hypothesis, this does not mean that resources are not important. The key message is that habitat changes resulting from decreasing light penetration are also important, and zooplankton consumers are responding to more than simply resources during long-term browning. These responses are also often taxon-or life history stagespecific, even for a single driver variable (e.g., compare across rows in Figure 7). Pronounced changes in the underwater light environment, which are particularly strong in clearwater lakes, alter vertical habitat gradients in UVR, PAR, temperature, and dissolved oxygen in ways that have important consequences for zooplankton physiology, habitat availability, and trophic interactions, including predation. Recovery from acidification and corresponding increases in pH, an integral part of browning in many lakes, may also interact with changes in the light environment in ways that alter not only total zooplankton abundance, but also their vertical distribution, trophic structure, and hence community structure. The limitations in predicting long-term trends from short-term data highlight the value of long-term data on the underwater light environment to accurately characterize changes in zooplankton communities in lakes undergoing browning.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study can be found in the Environmental Data Initiative (EDI, Williamson, 2019). The raw zooplankton data supporting the conclusions of this article are available upon request from the authors, and are in the process of being added to the EDI.

AUTHOR CONTRIBUTIONS
CW was involved in and supervised the long-term data collection and analysis, conceived of the study, and wrote the manuscript with input from, and sections written by EO, RP, and KW. EO, RP, and KW were also involved in the long-term data collection and analysis, with RP leading the statistical analysis. EO was the primary one to oversee quality control and archiving of data. All authors contributed to the manuscript revision, as well as read, and approved the submitted version.

FUNDING
Long-term data collection was supported by grants from the Andrew W. Mellon Foundation, Geraldine R. Dodge Foundation, and numerous National Science Foundation grants, including most recently NSF DEB LTREB 1754276. Data analysis, synthesis, and publication was supported by NSF DEB OPUS 1950170. We also thank the Miami University College of Arts and Sciences and the Ohio Eminent Scholar in Ecosystem Ecology fund.

ACKNOWLEDGMENTS
We thank Beth Mette and Eric Johnson for assistance in developing the database, Thomas Fisher for his statistical advice, and the students and colleagues, too numerous to mention by name, who assisted in the long-term data collection through the decades.