Annual winter water-level drawdowns in ﬂ uence physical habitat structure and macrophytes in Massachusetts, USA, lakes

. Annual wintertime water-level drawdowns are a common management strategy in recreational lakes; however, few studies have estimated their relative impact on lake littoral habitat among a set of typically co-occurring anthropogenic stressors including lakeshore development and herbicide application. Within 21 Massachusetts, USA lakes that represented a drawdown magnitude gradient (0.07 – 2.26 m), we assessed depth-speci ﬁ c littoral habitat (coarse wood, sediment, macrophytes) at two sites adjacent to forested or developed shorelines. Using generalized linear mixed models, we found coarse wood abundance and branching complexity was not correlated with drawdown magnitude but was primarily explained by the presence of lakeshore development. Drawdown magnitude was negatively correlated with silt cover and positively correlated with coarse substrate cover, with effects further varying by depth (0.5 m vs. 1 m). Macrophyte biomass and biovolume were negatively correlated with drawdown magnitude with effects also varying by depth for biomass. Macrophyte taxa with annual longevity strategies (e.g., Najas ﬂ exilis ) and amphibious growth forms increased in biomass proportions with drawdown magnitude. Distance-based redundancy analyses suggested macrophyte taxa composition was driven by drawdown magnitude, coarse substrate, alkalinity, water transparency, and herbicide use. These results suggest the importance of water quality and depth on macrophyte assemblage responses to winter drawdowns and the potential to develop drawdown-tolerant macrophyte assemblages. Altogether, understanding the unique impacts of anthropogenic stressors on littoral zone habitat across heterogeneous environmental lake conditions can help minimize impacts to lake ecological integrity while maintaining recreational value.


INTRODUCTION
Natural water-level fluctuations create spatiotemporal heterogeneity in the physicochemical habitat of lake littoral zones (Hofmann et al. 2008, Evtimova andDonohue 2015). Diverse littoral zone habitat (e.g., macrophytes, wood, bed texture) supports high within-lake diversity of invertebrates and fish (Weaver et al. 1997, Tolonen et al. 2001, White and Irvine 2003, provides fish spawning habitat (Winfield 2004, Lawson et al. 2011, mediates predator-prey interactions (Diehl 1992, Sass et al. 2006, Kornijów et al. 2016, contributes to whole-lake primary and secondary production (Vadeboncoeur et al. 2002, Vander Zanden et al. 2011, and may contribute to ecosystem resiliency (Kovalenko et al. 2012) by supporting longer food chains (Ziegler et al. 2015). In impounded systems, anthropogenic alterations to water level-alterations beyond the natural range of timing, magnitude, and frequency of daily to seasonal water-level fluctuations (Hofmann et al. 2008)-can impair the ecological integrity of littoral zones and hence lake ecosystems (Wantzen et al. 2008). Although scientific understanding of the role of natural (e.g., Evtimova and Donohue 2015) and modified (Leira andCantonati 2008, Zohary andOstrovsky 2011) water-level fluctuations in structuring littoral zone physical habitat has improved, there are limited empirical data on the impacts from prescribed water-level fluctuations regimes, including annual winter water-level reductions or drawdowns (referred to hereafter as winter drawdowns; Carmignani and Roy 2017).
Winter drawdowns are a widespread management practice conducted in impounded lakes in boreal and temperate regions typically as a consequence of power demands and flood protection in hydroelectric reservoirs (e.g., Mjelde et al. 2013), or as a strategy to reduce submerged macrophyte densities that may affect some recreational activities (Cooke et al. 2005). Drawdowns are initiated in fall and winter months, whereby water levels are reduced to desired minimum levels and rise to full pool levels upon spring flooding (Mattson et al. 2004). Through desiccation and accelerated erosional processes, drawdowns can reduce fine-textured sediment Matthews 2004, Cooley andFranzin 2008), organic matter, and nutrients , Furey et al. 2004 in exposure zones, leaving behind primarily larger sediment particles with low nutrient storage capacity. These abiotic changes along with direct physiological stresses from desiccation and freezing conditions can reduce macrophyte abundance and alter assemblage composition within drawdown exposure zones (Wilcox and Meeker 1991, Wagner and Falter 2002, Turner et al. 2005. Specifically, winter drawdowns can reduce macrophyte species reliant on vegetative structures for future propagation (i.e., perennials) in favor of high seed-bearing taxa (i.e., annuals) or taxa with multiple viable propagation strategies (reviewed in Carmignani and Roy 2017). Ultimately, these littoral habitat changes with drawdown can result in less complex physical habitat structure with negative implications for invertebrate and fish assemblages Meeker 1992, Meeker et al. 2017).
Where winter drawdowns occur, they are typically not the only disturbance contributing to loss in littoral zone habitat complexity (Kaufmann et al. 2014); lakeshore development, herbicide application, and nutrient loading also alter littoral habitat in drawdown lakes. Impacts from lakeshore development include reduced coarse wood (Christensen et al. 1996, Francis andSchindler 2006), reduced emergent and floatingleaved vegetation (Radomski and Goeman 2001, Alexander et al. 2008, Hicks and Frost 2011, finer sediments (Jennings et al. 2003), and lower sediment organic matter content (Francis et al. 2007). Lake nutrient enrichment in combination with other pressures that affect food web dynamics (e.g., fish winterkills, invasive species) can enable declines of submerged macrophytes particularly in shallow lakes (Phillips et al. 2016). However, disentangling the individual and potentially collinear effects of these anthropogenic stressors can be challenging (Van Sickle 2013), and elucidating the interacting effects of winter drawdowns with co-occurring anthropogenic stressors offers a novel area for research.
We aimed to determine the effects of winter drawdowns on physical habitat (i.e., coarse wood, sediment, macrophytes) of the littoral zone for lakes with decades of annual winter drawdowns. Given that littoral zone physical habitat can exhibit substantial inter-lake variability (Gasith andHoyer 1998, Weatherhead and, our study included 21 lakes that encompass a gradient of drawdown magnitude while attempting to account for other environmental gradients (e.g., water chemistry, morphometry) and co-occurring anthropogenic pressures (e.g., local shoreline development, herbicide application) that influence physical habitat. Our study can help refine adaptive lake management strategies to minimize ecological impacts in the context of multiple anthropogenic stressors.

Lake selection and study area
We selected impounded lakes (referred to hereafter as lakes) in the state of Massachusetts (MA), USA, using a stratified random approach to primarily capture a winter drawdown magnitude gradient. Lakes were selected from local conservation commissions and lake associations that responded to a statewide email survey (i.e., 397 out of 2080 waterbodies). We targeted lakes in the Northeastern Highlands (e.g., Western New England Marble Valleys/Berkshire Valley/ Housatonic and Hoosic Valleys) and two ecoregions in the Northeastern Coastal Zone (e.g., Connecticut River Valley, Lower Worcester Plateau) to help reduce water chemistry variation among waterbodies based on watershed land cover and geology ( Fig. 1; Griffith et al. 2009).
Where we received reported drawdown magnitude information (n = 21 lakes), we selected two lakes each from four drawdown magnitude classes (<0.5, 0.5-1, 1-1.5, >1.5 m) to ensure a drawdown magnitude gradient. We then selected eight additional lakes with a history of annual winter drawdowns but without magnitude information that were stratified into four lakeshore development density classes (e.g., 0-155, >155-284, >284-395, 412-536 buildings/km 2 ) calculated within a 100-m buffer around shore and determined by natural breaks in the data distribution. The final four lakes had no history of annual winter drawdowns, and these lakes were randomly selected based on lake area (0.012-0.073 or 0.11-0.89 km 2 ) and lakeshore development density (<97 or >105 buildings/ km 2 ). Where waterbodies were exhausted within a stratification (low drawdown magnitude class v www.esajournals.org <0.5 m), we extended our selection area to include the New England Coastal Plains and Hills in eastern MA, and randomly selected Silver Lake and Lake Boon (Fig. 1). We were unable to sample five of the original 20 selected lakes in 2014 due to access issues and replaced those with six additional lakes that are within our study area and represent lakes with current drawdown regimes or with no history of annual winter drawdowns, for a total of 21 lakes (Table 1).
Study lakes were in the Northeastern Highlands and Northeastern Coastal Zones (level 3 ecoregions) located in the Housatonic, Connecticut, Thames, Merrimack, and Blackstone River watersheds (Fig. 1). Inland MA has a continental temperate climate with four seasons. Mean minimum and maximum July and January temperatures for ecoregions in the Northeastern Highlands tend to be 1°-3°C degrees lower than in Northeastern Coastal Zone (Griffith et al. 2009). Winter precipitation averages 21. 6-25.4 cm (1981-2010) across the study area (National Oceanic and Atmospheric Administration, https:// www.ncdc.noaa.gov/cdo-web/datatools/norma ls). Lake watersheds have mixed land use with variable urban development ranging from 2% to 40% (median = 9%) with a general increase from west to east, and relatively small proportions of pasture (0-15%) and agriculture (0-8%). Total watershed forest cover ranged from 20% to 83% (median = 64%) among lakes. Forests are primarily composed of mixed deciduous and conifer stands including northern, central, and transition hardwoods. Lakes located in the Northeast Highlands are characterized by coarse-loamy to loamy soils and metamorphic Notes: Lake locations are mapped in Fig. 1. Abbreviations are non-WD, non-winter drawdown lake; WD, winter drawdown lake; NK, data is not known; -, data not applicable.
† Indicates non-drawdown lakes such that drawdown magnitude represents average low winter water levels.
Italic values indicate mean, minimum, and maximum values for drawdown magnitude, total phosphorous (TP) and lakeshore residential development.
v www.esajournals.org bedrock or limestone derived coarse-loamy soils and calcareous bedrock. In the Northeast Coastal Zone, lakes are underlain with sedimentary bedrock and alluvium soils, metamorphic bedrock with coarse-loamy soils, or coarse-loamy and sandy soils (see Griffith et al. 2009 for more detail).

Physical habitat sampling
We sampled lakes once in 2014 (n = 15 lakes) or 2016 (n = 6 lakes) in July-August when water levels were at or near full pool and macrophytes were generally at peak biomass. Since annual drawdown regimes have been maintained for at least two decades (Table 1), our single season sampling was presumed to reflect a sustained drawdown effect. At each lake, we established two sampling sites that stretched along 20-m shoreline segments. One site was selected with predominant forest riparian cover and the other site by human development (i.e., houses, lawns), each buffered by 50 m of similar shoreline land cover composition on each end. Sites were selected to represent shorelines sheltered from predominant wind-wave action and with gently graded slopes (i.e., ≤10%) to ensure we sampled conditions that support macrophyte biomass (Duarte and Kalff 1990).
We aimed to capture the major physical littoral habitat components including coarse wood, sediment, and macrophytes. At the site level, we enumerated all coarse wood (i.e., wood ≥10 cm in diameter at its thickest cross-section) at depths ≤1 m along 100 m of shoreline centered around the 20-m sites. Using methods from Newbrey et al. (2005), we quantified the branching complexity for each coarse wood piece. Branching complexity is an index based on the summation of inverted Strahler stream order (Strahler 1957) such that the bole or trunk has an order of 1, each extending branch from the bole an order of 2, and so on (see Newbrey et al. [2005] for detailed calculation). Complexity values can range hypothetically from 1 to infinity, but observed values were <400. For every site, we set three transects spaced 10 m apart and perpendicular to shore that extended to 1.5-2 m depths. Along each transect, we collected habitat data at 0.5, 1 m, and between 1.5 and 2 m depth contours classified as >1 m. Using a 1-m 2 quadrat, we visually estimated percentages of submerged macrophyte cover and biovolume, and sediment size classes (e.g., silt, sand, gravel, pebble, cobble, boulder). We defined macrophyte biovolume as the approximate percent of macrophytes that filled the water column within the sampling quadrat.
To determine macrophyte biovolume, we measured the approximate mean plant height, divided it by water depth, and multiplied by macrophyte cover. We summed the gravel, pebble, and cobble sediment size-class proportions per quadrat to create an aggregate coarse sediment variable to attain more non-zero data for analysis.
For sites sampled in 2014 (n = 15), we collected triplicate samples of the top 2 cm of sediment using 50-mL falcon tubes adjacent to a randomly selected 1-m 2 quadrat at each depth and site. Sediment samples were put on ice, kept frozen in the laboratory before percent organic matter content determination. Sediment was dried at 60°C for ≥24 h, weighed, placed in a loss-on-ignition furnace for 4 h, and weighed again to determine percent organic matter content. Depth-specific samples <1 g were aggregated.
Within the 1-m 2 quadrat, we randomly placed a 0.25-m 2 quadrat, harvested the above-ground portion of macrophytes within the smaller quadrat, and brought the macrophytes to the laboratory for identification and biomass measurement. Macrophytes were identified to species using Crow and Hellquist (2000a, b) except for Utricularia species and macroalgal taxa Chara and Nitella, which were left at genus. Individual macrophyte taxa were dried at 60°C for ≥24 h and weighed. Quadrat-level data were averaged across transects for each depth contour per site for statistical analysis.
We assigned macrophyte taxa to functional trait states based on morphology, longevity, amphibious capacity, fecundity, and native or non-native status (Appendix S1: Table S1). Previous studies have suggested these traits are influenced by annual winter drawdown regimes Meeker 1991, Cooke et al. 2005) and other water-level fluctuation disturbances (Willby et al. 2000, Arthaud et al. 2012. Taxa were assigned morphology states based on leaf arrangement and general plant height following nomenclature from Wilcox and Meeker (1991) and Meeker et al. (2017) including erect-and low-growth aquatics, low-rosette, and matv www.esajournals.org former. Low growth was defined as taxa that typically do not extend throughout the majority of the water column relative to other taxa (i.e., erect-aquatics). Longevity was categorized into perennial and annual taxa, along with perennials and annuals that possess storage organs (e.g., dormant buds in annuals, see Grime et al. 1990, Willby et al. 2000, Combroux et al. 2001, Hill et al. 2004, Capers et al. 2010, Arthaud et al. 2012). We divided taxa as amphibious or not following Willby et al. (2000); we expected amphibious taxa to be more tolerant of drawdown exposure. Lastly, fecundity was based on the number of reproductive organs (low < 10, medium = 10-100, high = 100-1000 yr −1 Áindividual −1 ) and divided by mode of reproduction as only seeds or as seeds and vegetative propagules following Willby et al. (2000) and Arthaud et al. (2012). We expected annuals and/or taxa with high reproductive output or multiple propagation strategies to be more tolerant of winter drawdowns. Species native status was determined using the PLANTS database (https://pla nts.sc.egov.usda.gov/java/index.jsp) and GoBotany databases (https://gobotany.nativepla nttrust.org/). If we could not locate trait information for taxa, we used descriptions from taxonomic keys (e.g., Crow and Hellquist 2000a, b; PLANTS database).

Water quality
We sampled water quality and determined secchi depth at the deepest part of each lake for two years between 2014 and 2017. In June, July, and/or August, we collected surface water samples for total phosphorous (TP), total nitrogen (TN), alkalinity, and dissolved organic carbon (DOC) that were analyzed at the University of New Hampshire Water Quality Analysis Laboratory. TP and TN were directly sampled with acid-washed polyethylene bottles, frozen, and analyzed through alkaline persulfate digestion followed by colorimetric measurement for PO 4 and NO 3 , respectively (Patton and Kryskalla 2003). Water samples for alkalinity and DOC were filtered through a pre-ashed microfiber glass filter, put on ice, cooled, and kept frozen, respectively. DOC was measured using US EPA (1979) with high temperature catalytic oxidation and alkalinity using the inflection point titration method.
Chlorophyll-a was filtered using a pre-combusted microfiber glass filter, put on ice, and kept frozen for <2 weeks before processing at the University of Massachusetts Amherst. We followed EPA method 445.0 in vitro determination of chlorophyll-a by fluorescence. Briefly, chlorophyll was extracted from the filters using 90% acetone with 18-24 h of extraction time. Extracted chlorophyll was measured using an AquaFluor fluorometer (Model 8000-010; Turner Designs, Sunnyvale, California, USA) and then acidified using hydrochloric acid to determine chlorophyll-b. Chlorophyll-b values were backcalculated to determine chlorophyll-a concentration in the original sample volume (Arar and Collins 1997).

Lakeshore development, herbicide use, and fetch
At the lake-level, we used the 2011-2014 Mass-GIS Building Structures (2-D) data layer to estimate lakeshore residential density as the number of buildings within a 100-m buffer around the shoreline. At the site level, we estimated effective fetch following methods from Håkanson and Jansson (1983) and Cyr et al. (2017). Over-water distances were measured in ArcGIS 10.3.1. Wind speeds and directions were taken from the United States National Oceanic and Atmospheric Administration using daily wind from Orange Municipal Airport, MA (USW00054756) running from 1998 to 2017. Our study lakes variably undergo herbicide application for nuisance macrophyte species during spring and summer seasons (Table 1). We assigned the presence or absence of herbicide application over the past two years for each site within each lake using annual herbicide use reported to the Massachusetts Department of Environmental Protection.

Lake water levels
We continuously monitored water levels for each lake from September/October of 2014 or 2015 to December 2017. We installed paired nonvented pressure transducers (Onset HOBO U20L-01, Bourne, Massachusetts, USA) at the point of outflow underwater and above water on shore and were both set to record at 2-h intervals. Paired pressure measurements were converted to water levels using HOBOWarePro (version 3.7.8, Onset Computer Corporation, Bourne, Massachusetts, USA). To calculate drawdown magnitude, we first isolated drawdown events using daily means by identifying the drawdown initiation date as the first record of consistent water level decline in the fall (i.e., October--November) and drawdown end date as the first record reaching pre-defined summer pool levels in winter-spring (i.e., drawdown end in January-June). We identified summer pool levels (i.e., drawdown refill target) as the median water level from non-drawdown phases in 2015 (n = 15) or as the top of spillways (i.e., lake water levels at full capacity without overflow downstream, n = 6). We determined drawdown magnitude as the lowest water level during drawdown relative to summer pool levels and used the average from the 2-3 drawdown events per lake for analyses.

Statistical analyses
We analyzed habitat response variables (macrophyte biomass, macrophyte biovolume, silt-sized sediment, coarse-sized sediment, percent organic matter, coarse wood abundance, coarse wood complexity) using generalized linear mixed models to fit various probability distributions and account for non-independence inherent in our nested study design (Appendix S1: Table S2, Bolker et al. 2009, Zuur et al. 2009). Macrophyte biomass did not fit a normal (Shapiro-Wilk, W = 0.41, P < 0.001) or log-normal error distribution (Shapiro-Wilk, W = 0.95, P < 0.001); hence, we used a gamma distribution with a log link and transformed the data using x + 0.001 g to elevate zero-values. We modeled percent sediment organic matter, macrophyte biovolume, and sediment size proportional data using a beta error distribution with a logit link, and applied the transformation derived from Smithson and Verkuilen (2006) to meet beta error distribution range values >0 and <1. We modeled site total coarse wood abundance and branching complexity count data by applying a negative binomial error distribution with a log link and an offset of coarse wood abundance for branching complexity counts.
We anticipated habitat responses to covary by sample depth (e.g., 0.5, 1, >1 m) along our drawdown magnitude gradient (Table 1), because of variable drawdown exposure and independent effects of depth on habitat. Thus, contour-level habitat response variables (i.e., all except coarse wood variables) were modeled with a drawdown magnitude-depth interaction, other potential environmental covariates, and lake as a random intercept (Appendix S1: Table S2). Since sediment organic matter was sampled in a subset of lakes (n = 15) and can potentially influence macrophytes, we also developed a separate set of models for macrophyte biomass and biovolume with organic matter as a predictor. We also applied generalized linear mixed models to each macrophyte trait state with sufficient non-zero values across the drawdown magnitude gradient using the same predictor structure as macrophyte biomass and biovolume models. Models were not applied to annuals with storage organs (longevity), moderate and high numbers of reproductive organs with seeds only (fecundity), mat-former and low-rosette (morphotype), and for non-native taxa (status). For coarse wood abundance and branching complexity, we tested an interaction between drawdown magnitude and lakeshore type (e.g., forested/developed).
We started with full predictor sets (Appendix S1: Table S2) of known covariates that could affect habitat response variables and iteratively removed single nonsignificant (P > 0.05) predictors using chi-square tests to simplify models and isolate important predictors. All continuous variables were Z-scored transformed before analyses. We checked for covariate collinearity using scatterplot matrices (e.g., Pearson r < 0.7) for continuous predictors, and generalized variance inflation factors (e.g., GVIF < 3) among continuous and categorical covariates using the car package in R (Fox and Weisberg 2011, version 2.1-5). We found secchi depth was strongly correlated with DOC (r = −0.76) and chlorophyll-a (r = −0.70), and consequently included only secchi depth in our models. We compared models using corrected Akaike information criterion (AIC c ) to determine the most parsimonious and plausible models for each habitat response variable (Burnham and Anderson 2004). We further determined marginal R 2 values for generalized linear mixed models (marginal R 2 GLMM ) to estimate proportion of variance explained by fixed predictors using the delta method to estimate observation level variance (Nakagawa et al. 2017). Marginal R 2 GLMM values were calculated with the performance package in R (Lüdecke et al. 2020, version 5.0). Models were validated by examination of residual plots at v www.esajournals.org predictor and model levels to ensure no patterns existed. We generated all regression models using the glmmTMB package (Brooks et al., 2017, version 0.2.1.0) performed in R (R Core Team, 2017, version 3.4.2).
We used distance-based redundancy analysis (db-RDA) to assess relationships between macrophyte taxa composition and winter drawdown magnitude and other environmental variables. db-RDA can use non-Euclidean distance measures more suitable for species community composition data (Legendre and Anderson 1999). We used Bray-Curtis to represent dissimilarity in species composition across sites (Bray and Curtis 1957), with the response matrix as macrophyte biomass at the contour level for 0.5 and 1 m depths. This yielded 84 samples (i.e., 21 lakes, 2 sites/lake, 2 contours/site). Before analysis, we first dropped rare taxa and traits with fewer than five observations (n = 20) and sites with no macrophyte biomass (n = 10), and subsequently logged transformed biomass of the remaining 21 taxa. To help isolate a set of environmental predictors to use in db-RDA, we used stepwise forward selection that minimizes AIC starting with the same full set of predictors in the univariate regression models. The final predictor set included drawdown magnitude, coarse substrate proportion, herbicide use, TP, alkalinity, and secchi depth. We also ran partial db-RDA to determine the relative influence of winter drawdown magnitude and coarse substrate on macrophyte taxa composition while controlling for other environmental variables (e.g., alkalinity, secchi depth). We performed an ANOVA-like permutational tests to assess whether environmental variables together (i.e., full ordination), each constraining axes, and individual environmental variables explained a significant proportion of variation in macrophyte taxa variation. The db-RDA constrained ordinations and permutational tests were conducted using the vegan package (Oksanen et al. 2019, version 2.5-3) in R.
Drawdown magnitude showed a marginally nonsignificant positive trend with wood complexity (β = 0.38, SE = 0.23, P = 0.099, Model R 2 GLMM = 0.46); however, this trend was driven by a forested site at the lake with the deepest drawdown (Otis) that had extremely high wood complexity. Wood was also less complex along developed shorelines than forested shorelines after accounting for coarse wood abundance (β = −0.87, SE = 0.39, P = 0.025; Table 2). We also found a positive effect of whole-lake residential density on wood complexity (β = 0.63, SE = 0.22, P < 0.001).

Sediment
Silt and coarse sediment proportions were moderately correlated with each other (Pearson r = −0.61), and this was reflected with similar predictor sets in our models (Table 2). Depth was significantly correlated with both silt and coarse substrate, whereby silt increased with depth and coarse particles decreased with depth. Silt proportion was best explained by an interaction between depth and drawdown magnitude (Table 3), whereby silt cover significantly decreased with drawdown magnitude at the 0.5 m depth (Fig. 2a). The top model for silt also included bed slope (steeper slopes had less silt), and lakeshore type (less silt in developed than forested sites) was included as a predictor in the next plausible model (Table 2). Coarse substrate was best predicted by the drawdown magnitude-depth interaction (Table 2), whereby coarse substrate significantly increased with magnitude at 0.5 and 1 m depths, with this effect waning with increased depth (Fig. 2b, Table 3). Organic matter content was significantly lower along developed shorelines and steeper slopes (Table 3). Drawdown magnitude showed nonsignificant negative effects on organic matter content (Fig. 2c), and this effect was strongest at the 0.5 and >1 m depth contours.

Macrophyte biomass and biovolume
Macrophyte biomass varied by 2-3 orders of magnitude, with mean biomass ranging from 0.17 to 73.44 g among lakes. The top model included a drawdown magnitude-depth interaction, alkalinity, secchi depth, and coarse substrate. Models with the addition of lakeshore type (developed/forested) and slope as predictors were also equally plausible models (i.e., <2 ΔAIC c ; Table 2). We found a negative correlation of drawdown magnitude on macrophyte biomass and the strength of this effect varied by depth (Fig. 3a). At the 1 m depth, drawdown magnitude showed a significant negative effect on biomass, while magnitude showed nonsignificant negative effects at 0.5 and >1 m depths (Table 3). Secchi depth and alkalinity had significant positive effects on macrophyte biomass, while coarse substrate was negatively correlated with macrophyte biomass (Table 3). The addition of organic matter as a predictor within a subset of lakes did not affect our interpretation on effects of winter drawdowns but had a significant negative effect on biomass (β = −0.52, SE = 0.23, P = 0.021).
Macrophyte biovolume varied from 1.1% to 34% among lakes. The top biovolume model was Notes: K represents the number of parameters and model weights are derived from models from full predictor sets to the top model. Abbreviations are RandI Lake , random intercept of lake; Mag, drawdown magnitude; Z, water depth; Alka, alkalinity; Csub, coarse substrate; ShoreType, lakeshore type (developed/forested); Herb, herbicide use (presence/absence); ResDens, shoreline residential density; Fetch, effective fetch; TP, total phosphorous; Secchi, secchi depth; OM, organic matter content; CWD, coarse wood abundance; Mag × Z, magnitude-depth interaction. similar to biomass (i.e., included a drawdown magnitude-depth interaction, alkalinity, secchi depth, and coarse substrate) with the addition of lakeshore type, whereby biovolume was lower along developed shorelines (9.1% AE 14) than forested shorelines (16% AE 16, Fig. 3b). Other plausible models included a negative effect of TP ( Table 2). As with macrophyte biomass, drawdown magnitude had a negative effect on macrophyte biovolume (Fig. 3b), which was significant at the 1 m depth and nonsignificant at 0.5 and >1 m depths (Table 3).

Macrophyte taxa and trait composition
Univariate response models for macrophyte traits showed variable responses to drawdown magnitude. For longevity traits, drawdown magnitude had no effect on perennials at 0.5 m depths but showed a marginally insignificant negative correlation at the 1 m depth (Appendix S1: Table S4). Also, the proportion of perennials were lower at 1 m compared with 0.5 m depths. In contrast to perennials, drawdown magnitude was positively correlated with annuals at 0.5 and 1 m depths, with a stronger effect at the 1 m RandI Lake 0.04 ---0.27 -0.30 -0.084 -Notes: Model terms include estimates (β) and standard errors (SE) for drawdown magnitude at 0.5, 1, and >1 m water depths, depth contrasts (e.g., 1-0.5 m), drawdown magnitude-depth slope contrasts (i.e., interactions), and other environmental covariates. Abbreviations are Secchi, secchi depth; Alka, alkalinity; Csub, coarse substrate; Slope, bed slope; ShoreType, contrast of developedforested shorelines; RandI Lake , random intercept of lake. Absence of a random lake intercept indicates a negligible variance term (e.g., <0.001). Bolded values indicate a significant correlation at P < 0.05. R 2 GLMM represents marginal R 2 values for fixed effects.
depth. Further, the proportion of annuals was higher at 1 m vs. 0.5 m depths and was positively correlated with alkalinity and herbicide use. Fecundity trait and morphotype proportions were not significantly correlated with drawdown magnitude or a drawdown magnitude-depth interaction (Appendix S1: Table S4). We found significantly lower proportions of the erect-growth aquatic morphotype at the 1 m depth compared with the 0.5 m depth and found the converse for low-growth aquatics. The proportion of amphibious taxa was positively correlated with drawdown magnitude at the 0.5 m depth (Appendix S1: Table S4). Additionally, we found higher amphibious proportions at 0.5 m compared with 1 m depths, and with higher v www.esajournals.org alkalinity, higher effective fetch, with less coarse substrate, and the absence of herbicide use (Appendix S1: Table S4).
The partial db-RDA of drawdown magnitude and coarse substrate explained 9.0% of macrophyte composition, whereby drawdown magnitude was negatively correlated with db-RDA1 and coarse substrate negatively correlated with both axes (Fig. 4b). The partial db-RDA triplot helped to isolate apparent trends between macrophyte species and drawdown magnitude and coarse substrate while controlling for other important environmental gradients. For example, sites with relatively small drawdown magnitudes possessed high biomass of N. odorata, Utricularia species, Potamogeton robbinsii, and V. americana. In contrast, high biomass of Najas flexilis and Najas minor, Chara, Gratiola aurea, and Isoetes species corresponded to deeper drawdown magnitudes and coarser substrates (with the exception of Chara associated with finer substrates). Overall, lower biomass for the majority of species was found at sites with high proportions of coarse substrates.

DISCUSSION
We provide evidence that annual winter drawdowns alter littoral zone physical habitat even at relatively small magnitudes of <2 m. At depths typically exposed during drawdowns (e.g., 0.5 Macrophyte taxa scores are represented as abbreviated taxa codes (see Appendix S1: Table S1). Only environmental variables with inter-set correlation coefficients > |0.5| on the first two axes are shown. Abbreviated environmental vectors are (bolded) Csub, coarse substrate proportion; Mag, drawdown magnitude; Alka, alkalinity; Secchi, secchi depth. v www.esajournals.org and 1 m, depending on the lake), we found significant changes in sediment texture, macrophyte abundance, and macrophyte taxonomic and functional composition as a function of drawdown magnitude. Concordantly, at depths rarely exposed by drawdowns (i.e., ≥1.5 m), magnitude was not correlated with physical habitat components, suggesting that impacts from winter drawdowns correspond with the depth of exposure. Drawdown magnitude poorly explained coarse wood abundance and branching complexity variability; instead, coarse wood abundance and complexity was greatly reduced at developed shorelines compared with forested shorelines, demonstrating distinct effects of different anthropogenic activities on littoral zone habitat.

Winter drawdown effects on littoral habitat
Winter drawdowns coarsened sediment with associated reductions in silt cover and organic matter content at depths within exposure zones. These patterns are consistent with previous winter drawdown studies Falter 2002, Cooley andFranzin 2008) and other water-level fluctuation regimes (Evtimova and Donohue 2015) that suggest accelerated sediment focusing from exposure zones to depths below water-level minimums. As water levels decline, fine sediments at depths typically protected from wave action at normal water levels become susceptible to resuspension and are transported to deeper depths (Effler et al. 1998, Dirnberger andWeinberger 2005). Furthermore, water column mixing likely temporally overlaps with water-level recession from drawdowns in October to December, which may enhance sediment focusing (Effler and Matthews 2004). Ultimately, the likely interaction between annual drawdowns conducted for several decades and short-term high wind/ wave events (Hofmann et al. 2008) has coarsened exposure zones (Hall et al. 1999, Furey et al. 2004. We found annual winter drawdowns affect the abundance, taxonomic, and functional composition of submerged macrophytes in drawdown exposure zones. Consistent with previous winter drawdown studies (Siver et al. 1986, Turner et al. 2005, Olson et al. 2012, measures of macrophyte abundance (e.g., biomass and biovolume) were negatively correlated with drawdown magnitude, particularly at the 1 m depth. Drawdowns did not affect macrophyte abundance at depths >1 m, presumably because they are rarely exposed during drawdown, and at the 0.5 m depth because other environmental factors (e.g., ice erosion, Renman 1989, Hellsten 1997, may be more important at shallow depths. The correlations between drawdown magnitude, coarse substrate, and macrophyte biomass suggest winter drawdowns reduce macrophytes directly through exposure to winter conditions and indirectly through sediment coarsening over time. Wagner and Falter (2002) similarly found significantly lower macrophyte biomass on cobble substrate, which existed at higher frequencies in shallow-exposed depths in an annual winter drawdown lake. Macrophyte abundance tends to decrease with increasing sediment particle size (Anderson and Kalff 1988) because of low nutrient diffusion rates and nutrient capacity (Barko and Smart 1986), and its association with relatively high wind/wave energy and steeper littoral slopes Kalff 1986, Cyr 1998). Furthermore, winter drawdowns may decouple positive feedbacks between macrophyte beds, fine sediment accretion, and erosional reduction (Barko and James 1998), and enable sediment coarsening and further macrophyte reduction over time.
Taxa that appeared to be sensitive to winter drawdowns were Nymphaea odorata, Vallisneria americana, and Potamogeton robbinsii. Previous studies have also shown declines of P. robbinsii (Beard 1973, Nichols 1975, Crosson 1990) associated with winter drawdowns. These species are perennial taxa that primarily propagate via vegetative structures (e.g., rhizomes), which have been hypothesized to be sensitive to desiccation, freezing, and erosional disturbance related to winter drawdown (Rørslett 1989, Wagner andFalter 2002). Accordingly, we found a decline in perennial taxa, particularly at the 1 m depth including the aforementioned species. We found proportionally more perennials at the 0.5 m depth compared with 1 m depth and no effect of drawdowns at 0.5 m, suggesting that perennial taxa are variably susceptible to winter drawdown disturbance. For example, we found relatively high proportions of perennial taxa that are also amphibious (Gratiola aurea and Sagitarria) at the 0.5 m depth in high drawdown magnitude lakes and were generally less abundant at the v www.esajournals.org 1 m depth. Perennial taxa have plastic and variable propagation strategies (Barrat-Segretain et al. 1998, Combroux andBornette 2004), high niche breadth (Alahuhta et al. 2017), and ability to colonize exposure zones late in the growing season (August-September). Furthermore, the inter-annual variability of drawdown exposure weather conditions (e.g., freezing temperatures, snowfall) could permit variable rhizome survival (Lonergan et al. 2014).
Winter drawdowns can select for drawdowntolerant macrophyte assemblages (Siver et al. 1986, Richardson et al. 2002. Where macrophytes were present, several taxa were positively associated with drawdown magnitude. Consistent with other studies, we found positive associations of N. flexilis (Beard 1973, Nichols 1975, Tazik et al. 1982, Crosson 1990, Turner et al. 2005 and N. minor (Siver et al. 1986), and the macroalgae Chara (Wagner and Falter 2002) with drawdown magnitude. These taxa generally possess an annual longevity strategy that are largely dependent on sexual diaspores in the form of seeds (Najas species) or oospores (Chara). Concordantly, drawdown magnitude was positively related to annuals at exposed depths, consistent with ruderal life history strategies (Grime 1977, Rørslett 1989. We also found a positive, albeit weak correlation between amphibious taxa (G. aurea, Sagitarria, Elatine minima) and drawdown magnitude at the 0.5 m depth, aligning with previous work (Rørslett 1989), although effects may be stronger under deeper drawdown magnitudes.
Several macrophyte traits were unrelated to drawdown magnitude. We observed no correlation between drawdown magnitude and taxa with moderate to high fecundity levels that produce both seeds and vegetative propagules, a finding consistent with Arthaud et al. (2012), suggesting several reproductive strategies may enable a taxa's persistence in annual drawdown regimes. We also found no distinct trends among macrophyte morphologies and drawdown magnitude. Previous studies found increases in matforming and low-rosette taxa with drawdowns (Wilcox and Meeker 1991); however, our dataset was insufficient to assess changes in these morphologies because of low sample sizes. Wilcox and Meeker (1991) also found declines in lowand erect-growth aquatics with drawdowns; the lack of a relationship in our study may be explained by our relatively mild amplitudes in combination with co-occurring alkalinity and secchi gradients.
Effects of lakeshore development, water quality, and herbicide use on littoral habitat In addition to winter drawdown magnitude, lakeshore development contributed to explaining littoral habitat variability. Along lakeshore development, our observed patterns of reduced and less complex coarse wood (Christensen et al. 1996, Jennings et al. 2003, Francis and Schindler 2006, Merrell et al. 2009), reduced sediment organic matter (Francis et al. 2007), and reduced macrophyte abundance and altered macrophyte composition (Radomski and Goeman 2001, Jennings et al. 2003, Merrell et al. 2009 Vondracek 2017) has been well documented. Decreases of coarse wood recruitment is a function of lake riparian deforestation and direct removal of coarse wood within a lake (Francis and Schindler 2006). Lower structural complexity along developed shorelines compared with forested sites may be due recreational driven processes such as wave erosion from motorboats, physical removal of branches for firewood, or to reduce angling interference (Newbrey et al. 2005). In contrast, whole-lake residential density was positively correlated with branching complexity after accounting for coarse wood abundances. For several lakes with relatively high lake-wide residential development, high branching complexity resulted from many small wood pieces accumulated around coarse wood, particularly along forested shorelines. Often, forested shorelines in these lakes exist only along windsheltered coves and/or adjacent to inlets where coarse and smaller pieces of wood may collect over time and where beaver populations may colonize (Marburg et al. 2006). Deforested lakeshores and less coarse wood via lakeshore development can lessen organic matter retention particularly at shallower depths (Francis et al. 2007, Merrell et al. 2009). Consequently, existing organic matter may be transported to deeper depths via erosional forces from wave action and drawdown, which matches reported depth distributions associated with lakeshore development (Francis et al. 2007). As with coarse wood, macrophytes are directly removed (Asplund andCook 1997, Radomski andGoeman 2001) via management strategies (e.g., hand-pulling, herbicide, suction and mechanical harvesting, benthic barriers) to facilitate recreational activities, particularly in front of active lakefront property (Payton and Fulton 2004). Submergent taxa were the dominant growth form and we detected emergent and floating-leaf taxa in only 5% and 17% of our sampling quadrats, as seen in previous studies (Radomski and Goeman 2001, Jennings et al. 2003, Dustin and Vondracek 2017. Despite the dominance of submerged taxa, we found lower macrophyte biovolume along developed vs. forested shorelines, supporting previous observations of macrophyte cover (Merrell et al. 2009). This effect may correspond to less tallgrowing submerged taxa in our study that may impede recreation.
Water quality factors also influenced macrophyte composition and total abundance metrics. Macrophyte biomass and biovolume were positively correlated with alkalinity. The biomassalkalinity trend supports previous observations (Duarte and Kalff 1990) and the positive correlation between biovolume and alkalinity may result from relatively low-lying species within the water column (e.g., isoetids) associated with low alkaline lakes along with higher biomass in more alkaline lakes. Alkalinity is a major environmental factor controlling macrophyte species composition (Roberts et al. 1985, Vestergaard and Sand-Jensen 2000a, Alexander et al. 2008) because of its tight correlation with bicarbonate (HCO 3 − ) concentrations that can be variably used as a carbon source for different macrophyte species (Madsen and Sand-Jensen 1991). Higher alkaline lakes tend to support more macrophyte species (Roberts et al. 1985) composed predominantly of the more species-rich elodeids and charophytes compared with soft-water lakes with more isoetids (Vestergaard and Sand-Jensen 2000b). We observed Chara, P. pusillus, Vallisneria americana, and Myriophyllum spicatum associated with moderate to high alkaline conditions and Nitella, N. odorata, B. schreberi, Isoetes, Utricularia, and Potamogeton bicupulatus associated with low alkalinities, which is consistent with previous studies (Alexander et al. 2008, Capers et al. 2010. Further, annual taxa were positively related to alkalinity, which likely derives from increased abundances of Chara beds in more alkaline conditions. Water transparency directly influences the amount of colonizable area for macrophytes where increases in clarity allow for deeper macrophyte colonization (Chambers and Kalff 1985, Duarte and Kalff 1990, Vestergaard and Sand-Jensen 2000b and increases in macrophyte biomass and cover (Barko et al. 1982, Cheruvelil andSoranno 2008). Low-lying species can persist at deeper depths in high clarity conditions (e.g., Isoestes, Mjelde et al. 2013), as we found for proportions of low-growth aquatic taxa. Although the effect of water clarity on abundance is typically more important at deeper depths (>2 m, Duarte and Kalff 1990), we were able to detect an effect because several lakes exhibited relatively low clarity (e.g., <2 m visibility). In our study, secchi depth was negatively correlated with DOC and chlorophyll-a, which influence water transparency (Canfield andHodgson 1983, Brezonik et al. 2019). Although the importance of specific drivers of water clarity variability is lakespecific, high chlorophyll-a concentrations (Kissoon et al. 2013) or DOC (McElarney et al. 2010 can limit depth range distributions and growth of submerged macrophytes. Herbicide use also structured macrophyte taxa composition. The lack of herbicide use overlapped with relatively higher alkalinity and secchi depths that together shaped a distinct macrophyte assemblage including the non-native invasive species M. spicatum. Interestingly, annual taxa were positively correlated with herbicide use. Annual taxa emerging from seed banks may become relatively abundant in the following growing season after targeted taxa are treated (Hussner et al. 2017). Although herbicide use influenced macrophyte composition, it was not an important predictor of total macrophyte biomass and biovolume compared with other environmental predictors including drawdown magnitude. Furthermore, the likely inter-lake variability of the type of herbicide, the associated selectivity to a target species, and the spatial pattern of application likely muted herbicidemacrophyte biomass patterns. Importantly, our study lakes with or without herbicide application were evenly distributed across the sampled drawdown magnitude gradient and did not confound our interpretation of drawdown-macrophyte patterns.

Implications for winter drawdown management
A primary reason for the implementation of annual winter drawdowns is to reduce nuisance densities of aquatic vegetation that inhibit recreational activities (Cooke et al. 2005). Our results show that drawdowns can partially meet this objective, as we observed a general decrease in macrophyte biomass and biovolume at depths exposed during drawdown across various ambient water quality conditions. However, macrophytes are not completely lost from exposure zones and considerable variability exists among lakes. Macrophytes can recolonize into exposure zones after a drawdown via seed banks or vegetative propagules from macrophytes at deeper unexposed depths and eventually resulting in a drawdown-tolerant macrophyte assemblage (e.g., Turner et al. 2005). Species that can rapidly colonize exposure zones upon refill are at an advantage over slow-growing species and can include potentially invasive species (Crosson 1990). The widespread invasive Eurasian milfoil (M. spicatum) is a frequent target of winter drawdowns, and we found relatively low biomass of M. spicatum in drawdown-exposed areas of 4 lakes, consistent with previous studies (Lonergan et al. 2014). This suggests drawdown can limit but not eliminate this species probably because of specific freezing and/or drying threshold conditions needed to prevent regrowth (Lonergan et al. 2014) and the ease of dispersal via fragmentation from unimpacted, deeper depths. Other invasive species tolerant to drawdown conditions, such as N. minor, may proliferate in drawdown exposure zones. After declines of M. spicatum from two winter drawdowns, Siver et al. (1986) observed increases in N. minor and N. flexilis in exposure zones in a Connecticut lake. Often, other macrophyte management strategies (e.g., herbicide application) are needed to supplement winter drawdowns to sufficiently control or eradicate target species over longer time periods (Cooke et al. 2005).
Our data suggest macrophyte responses to drawdown magnitude are likely modified by the environmental context in littoral zones and lakes. Winter drawdown regimes may impact macrophytes relatively more in littoral zones with low water clarity or low alkalinity than under high alkaline or high water clarity conditions where macrophyte colonization and biomass production can be extensive. Furthermore, lakes with high clarity or alkalinity may have a higher probability to develop a drawdown-tolerant macrophyte assemblage because of a richer species pool (Vestergaard and Sand-Jensen 2000b). Therefore, applying an equal drawdown magnitude across lakes with varying water quality conditions will have varying macrophyte impacts. Identification of winter drawdown-tolerant and sensitive taxa associated with different water quality conditions will require macrophyte surveys across many lakes within lake water quality classifications as seen in Mjelde et al. (2013) with oligotrophic and low alkaline lakes.

Conclusion
Multiple anthropogenic stressors degrade littoral zone habitat structure important for littoral zone biota (Miranda et al. 2010). In recreational lakes of Massachusetts, annual winter waterlevel regimes, lakeshore development, and herbicide application impact physical habitat through changes in littoral zone sediments, macrophyte assemblages, and coarse wood. Drawdown impacts are depth-specific and observed even at relatively mild drawdown magnitudes. Additionally, the variable state of macrophyte assemblages (i.e., tolerant taxa) in exposure zones suggests the importance of environmental context (e.g., water quality, spatial dynamics) at lake and watershed levels (e.g., land use) as seen in larger studies (e.g., Sass et al. 2010). Incorporating lake-specific, ambient environmental conditions into winter drawdown management can help to improve implementation of winter drawdowns while conserving ecological integrity. The alteration and reduction of complex littoral habitat can modify predator-prey interactions (Diehl 1992, Sass et al. 2006, Kornijów et al. 2016) and shape nutrient and energy flow in lake food webs (Barko and James 1998). Climate change will likely further affect littoral zone habitat availability through changes in lake water-level fluctuations. Summer drought conditions may become more frequent with climate change in the northeastern USA (Hayhoe et al. 2007) causing reductions in lake water levels and altering fish population dynamics (i.e., decreased fish growth) because of inaccessibility to critical spawning, predator refuge, and feeding habitat in littoral zones (Gaeta et al. 2014, Hardie andChilcott 2016). Limiting further habitat loss by protecting areas of complex habitat structure (e.g., inlets, forested shorelines) in these impaired lake ecosystems will be essential to preserve current ecosystem resilience to anticipated effects of climate change on lake water levels.