Watershed, lake and food web factors influence diazotrophic cyanobacteria in mountain lakes

Cyanobacterial blooms can occur in freshwater ecosystems largely isolated from development and not experiencing extensive cultural eutrophication. For example, remote mountain lakes can experience intense blooms of diazotrophic (nitrogen-fixing) cyanobacteria caused by factors acting at different spatial and temporal scales. In this study, we examined how cross-scale interactions among watershed, lake, and food web characteristics influence diazotrophic cyanobacteria biovolume in mountain lakes. We quantified diazotrophic cyanobacteria biovolume, zooplankton abundance, and physico-chemical variables for 29 lakes in the Cascade Mountains of Oregon, USA, in summer 2019. Watershed characteristics were compiled from historical data sets available for the region. Diazotrophic cyanobacteria biovolume ranged across the lakes from 0 to 1,930,000 µm3 mL-1; Dolichospermum was the most common genus. Random forest models showed that 11 watershed, lake, and food web characteristics explained 76% of the variance in diazotrophic cyanobacteria biovolume among the sampled lakes. Structural equation models suggested that the drainage ratio (i.e., relative area of the lake to the watershed) correlated positively with phosphorus concentrations and in turn, with diazotroph biovolume. Among lakes, hypolimnetic dissolved oxygen was negatively correlated with diazotroph biovolume, possibly due to the release of nutrients, like phosphate and iron, bound to sediments. In addition, zooplankton grazers were negatively related to diazotrophic cyanobacteria biovolume, potentially reflecting the influence of stocked fish.  Thus, lake management must account for bottom-up factors, such as nutrient loading, which is influenced by lake morphometry and watershed size, as well as top-down factors, such as fish stocking, to effectively mitigate diazotrophic cyanobacterial blooms.

Methods
We compiled a data set of 132 Oregon Cascade lakes on fish stocking, water chemistry, and watershed characteristics from various sources to identify a potential set of study lakes (D. M. Johnson 1985; USDA Forest Service 1996; US Environmental Protection Agency 2009; 2016).  We only considered lakes >10 hectares and with a maximum depth of ≥3 meters to exclude extremely small lakes or ponds that can be seasonally variable in depth and emergent vegetation cover (US Environmental Protection Agency 2009). We created a representative sample of lakes from this dataset using binary regression trees from the R package rpart v.4 (T. Therneau, Atkinson, and Ripley 2019). A binary regression tree repeatedly divides the response data into nodes to reduce variation within nodes based on predictor variables. Total P concentration was the response variable as P is often a crucial lake characteristic for diazotrophic cyanobacteria and varies significantly in the Cascades (D. M. Johnson 1985; D. P. Morris and Lewis 1988; J. J. Williams et al. 2016). In addition, total P was the most commonly available chemical analyte for all lakes from the compiled datasets.  We chose elevation, 10-year annual fish stocking average rate, and surface area as the predictor variables. From the resulting binary regression tree, we selected a similar number of lakes from each node to make up the 29 lakes of the study, which were sampled in summer 2019 (Appendix S1: Fig. S1). 

We compiled lake morphometry data from the Atlas of Oregon Lakes database, including bathymetry to calculate Schmidt Stability Index (a measure of stratification strength), using the R package rLakeAnalyzer v. 11 (D. M. Johnson 1985; L. Winslow et al. 2019). We obtained fish stocking records from 1978-2019 for the study lakes, including biomass, count, and species (Oregon Department of Fish and Wildlife, unpublished). The fish were mostly stocked as fingerlings that were < 250 mm long (fork length) and thus were primarily planktivores (D. A. Beauchamp 1990; J. J. Elser et al. 1995).

For many watershed variables, we extracted data from the LakeCat dataset, which provides watershed characteristics of lakes within the conterminous US derived from existing spatial data such as land cover, geology, and long-term climate (Appendix S1, Hill et al., 2018). Average slope and aspect of watersheds was derived from Digital Elevation Models with 10-m resolution (US Geological Survey 2019). Maximum snow-water equivalent (SWE) for each watershed in the previous winter (2018-2019) was derived from the Snow Data Assimilation System daily estimated SWE, selecting the days based on the maximum SWE at the nearest SNOTEL site for each watershed (National Operational Hydrologic Remote Sensing Center 2004; USDA Natural Resources Conservation Service 2020). Percent change in forest cover in each watershed for the past 20 years was derived from the spatial Global Forest Change dataset for 2000-2019 (Hansen et al. 2013). We calculated drainage ratio for each lake by dividing the watershed area by the lake surface area.

Lake & food web: field sampling

 Each lake (n=29) was sampled twice in summer 2019 to capture within season variation, with the first sampling occurring between June 24-July 13 and the second between August 8-31. We measured physical and chemical variables at the deepest spot in the lake, including a complete depth profile of temperature and dissolved oxygen with a YSI ProODO meter (Yellow Springs, Ohio, USA).  For the stratified lakes, the profile data were divided into the thermal layers (epilimnion, metalimnion, and hypolimnion) to calculate average epilimnion temperature, mixed layer depth (i.e. the bottom of the epilimnion), and average hypolimnetic dissolved oxygen. For unstratified lakes, average whole water column temperature was used in lieu of epilimnetic temperature and the bottom dissolved oxygen value was used in lieu of average hypolimnetic dissolved oxygen. 

We sampled for phytoplankton with a Van Dorn sampler (Wildco, Yulee, Florida, USA) at one meter below the surface at the deepest point in the lake, identified using bathymetric maps and a depth sounder. Samples were collected in 250-mL brown Nalgene bottles and preserved for identification using Lugol’s solution.  An additional sample was taken at the deep spot and a known volume was filtered onto glass fiber filters (Whatman GF/C, 1.2-m pore size) for chlorophyll-a analysis. We collected samples for nutrient analyses (total N and P) from the top 5 m of the water column, using a 5-m long, 2.54-cm diameter tygon tube and transferred them into 125-mL HDPE bottles, which were kept cool and preserved with H2SO4 until frozen in the laboratory. We conducted an integrated tow of the water column for crustacean zooplankton using a plankton net with 80-µm mesh and 25-cm diameter, starting two meters from the bottom and preserving the sample in 70% ethanol for identification.

Lake & food web: Lab analyses

 We concentrated preserved phytoplankton samples before counting by gently mixing the sample for 5 minutes and then taking a 100-mL subsample for settling in a graduated cylinder. After 100 hours of settling, the top 98 mL were removed via a vacuum pump and reserved to dilute while the remaining 2 mL were used for counting. We counted and identified 300 natural units per concentrated sample in Palmer counting cells to the genus level or to the lowest taxonomic level possible, using taxonomic guides (G. M. Smith 1950; J. D. Wehr 2002; R. A. Matthews 2016) with a Leica DM1000 microscope at 400X and ICC50 HD camera (Leica Microsystems Inc., Buffalo Grove, IL). Diazotrophic genera were determined based on current literature (Bergman et al. 1997; I. Berman-Frank, Lundgren, and Falkowski 2003; Reynolds 2006a). We measured the dimensions of 20 individuals of each taxon in each sample to calculate biovolume using standardized equations based on the shapes of taxa (Appendix S1: Table S4; Hillebrand et al., 1999). For zooplankton, we counted and identified 250 individuals from each preserved sample to the order level for Copepoda and to the family level for Cladocera using taxonomic guides (M. D. Balcer, Korda, and Dodson 1984; J. H. Thorp and Covich 2014) with a Leica M165C microscope at 100X and IC80HD camera. For further analyses, cladocerans were aggregated into two groups based on different feeding impacts, as Daphnia, large and efficient, and small cladoceran grazers, smaller and less efficient (e.g. Bosmina, Holopedium, Ceriodaphnia) (William R. DeMott 1982; Reynolds 2006a).

We extracted chlorophyll-a from filters using acetone for 20 hours in a dark refrigerator, and measured concentrations using a fluorometer following Arar & Collins (1997). We used a persulfate solution to digest total P samples heated to 100°C and then analyzed with a Shimadzu UV-1800 spectrophotometer (Kyoto, Japan) using the molybdenum blue colorimetric method (detection limit: 0.002 mg/L; precision limit: +/- 0.004 mg/L) (C. J. Patton and Kryskalla 2003; APHA 2018a). We also used a persulfate solution to digest total N samples heated to 100°C and then analyzed with a SmartChem 200 discrete analyzer (Guidonia, Italy) for colorimetric determination of nitrate and nitrite (detection limit: 0.01 mg/L; precision limit: +/- 0.01 mg/L) (APHA 2018b).