Bottom-Up Forces Drive Increases in the Abundance of Large Daphnids in Four Small Lakes Stocked with Rainbow Trout (Oncorhynchus mykiss), Interior British Columbia, Canada

The introduction of salmonids into lakes of western North America for sport fishing is a widespread phenomenon. While numerous investigations have documented cascading trophic interactions upon the introduction of fish into naturally fishless systems, little research has been done to investigate the importance of natural fish status (fishless vs. fish bearing) in modulating historical food web response to dual forcing by bottom-up (resource regulation from nutrients) and top-down (planktivory from stocked fish) processes. We used the paleolimnological record to reconstruct food web changes in four lakes in interior British Columbia that have been stocked with rainbow trout since the early to mid-1900s. Analysis of pigments, diatoms, and Cladocera was undertaken in cores from all lakes. We predicted that if fish were important in structuring cladoceran abundance and composition, we would document a decline in the abundance of large daphnids post-stocking in our two naturally fishless lakes, and little change in the two fish-bearing lakes. Instead, we documented increased abundance of large daphnids after stocking in all lakes in the early to mid-1900s, a finding inconsistent with size-selective predation from planktivorous fish. Further, our data suggest that deep, low-oxygen refugia may be important in sustaining populations of large Daphnia, a process which was enhanced by increased nutrients and lake production according to sub-fossil diatom and pigment analyses. This study shows that fish stocking does not invariably result in a decrease in large-bodied Cladocera and that nutrients and lake type can modulate the response of invertebrate planktivores.

tivory from stocked fish) processes. We used the paleolimnological record to reconstruct food web changes in four lakes in interior British Columbia that have been stocked with rainbow trout since the early to mid-1900s. Analysis of pigments, diatoms, and Cladocera was undertaken in cores from all lakes. We predicted that if fish were important in structuring cladoceran abundance and composition, we would document a decline in the abundance of large daphnids post-stocking in our two naturally fishless lakes, and little change in the two fish-bearing lakes. Instead, we documented increased abundance of large daphnids after stocking in all lakes in the early to mid-1900s, a finding inconsistent with size-selective predation from planktivorous fish. Further, our data suggest that deep, low-oxygen refugia may be important in sustaining populations of large Daphnia, a process which was enhanced by increased nutrients and lake production according to sub-fossil diatom and pigment analyses. This study shows that fish stocking does not invariably result in a decrease in

INTRODUCTION
The introduction of salmonids to lakes of western North America for sport fishing opportunity is widespread, with the earliest stocking events occurring in the mid-to late 1800s in western wilderness lakes. By 1988, 20-75% of mountain lakes in the western USA were stocked, many of which were naturally fishless (Bahls 1992). The most commonly introduced salmonids were rainbow trout (Oncorhynchus mykiss), brown trout (Salmo trutta), and brook trout (Salvelinus fontinalis) (Fuller and others 1999). The introduction of game fish can fundamentally alter food web structure, via both bottom-up (that is, resource regulation) and top-down forces (that is, predator regulation). For example, establishment of planktivorous fish populations have been associated with declines in zooplankton composition and size structure due to zooplanktivory, often characterized by declines in large-bodied (> 1 mm) Daphnia spp. (for example, Carpenter and others 1985;McNaught and others 1999;Parker and others 2001;Tiberti and others 2014). The loss of large and effective grazers can alter phytoplankton production by changing species composition (McNaught and others 1999) and increasing algal biomass (Sarnelle and Knapp 2005). Importantly, fish have also been suggested to directly alter nutrient cycling by excreting phosphorus (P) previously unavailable to planktonic communities (Leavitt and others 1994;Vanni and others 1997;Schindler and others 2001; Stuparyk and others 2019), although the relative importance of the mechanisms leading to increased phytoplankton production often associated with fish stocking is not always clear.
British Columbia (B.C.) maintains an active fish stocking program with over 800 lakes stocked annually, with the Thompson-Nicola region in the southern interior (SOI) experiencing the second highest angling activity in the province (FFSBC 2013). The lake types situated in the SOI differ from those of high-elevation mountain regions that have been investigated previously, in which many are meso-eutrophic (Ashley and Nordin 1999) and some contained native non-game fish prior to stocking due to easier colonization routes (fewer migratory obstacles). Additionally, because lakes within the B.C. interior are less remote than alpine basins, many sites have been subject to anthropogenic disturbances of the lakes and their watersheds, including impoundment, forest harvesting, and shoreline development (for example, Reavie and others 2000), all of which may favor nutrient delivery to lakes. Because nutrient loading may increase the amount of high-quality autochthonous carbon sources available to zooplankton (for example, Carpenter and others 2015) and/or alter phytoplankton species composition, shifts in both zooplankton biomass (for example, Yuan and Pollard 2018) and functional diversity related to resource acquisition (Barnett and Beisner 2007) may occur due to these bottom-up forces. Moreover, limnological changes often associated with eutrophication, such as the development of deepwater anoxia and a trend toward the predominance of cyanobacteria, can have important impacts on trophic interactions (Benndorf 1987;DeMott and others 2001;Hembre and Megard 2003).
Recently, some studies have focused on the impacts of introduced salmonids on lake biota in more productive lakes from lower-elevation regions of western North America. These studies have shown that naturally fish-bearing systems may exhibit little change in invertebrate species abundance, and composition and size structure following salmonid introduction (Nasmith and others 2012;Hanisch and others 2013). For example, microcrustacean zooplankton communities in naturally fish-bearing and stocked fish-bearing lakes have been found to be similar in size structure, abundance, and biomass, whereas fishless lakes are dominated by fewer but larger zooplankton (Holmes and others 2017). This pattern has been noted in other systems (Jeppesen and others 2000) and may be attributable to diet overlap between many native fish and stocked fish (Hanisch 2016). Although these findings have important implications for predicting how aquatic ecosystems will respond to stocking, less work has been done to investigate both top-down and bottom-up controls of variability in food web response to stocking of vertebrate planktivores in lakes of contrasting status (fish bearing vs fishless), especially over multiple decades. Given that lake sediments integrate primary and secondary production across all lake habitats, paleolimnological analysis of fossil zooplankton and phototrophic pigments can provide a holistic view of lake responses to fertilization and trophic manipulations over decadal-to-centennial timescales not amenable to modern experimentation (for example, Leavitt and others 1989).
The goal of this study was to investigate the importance of fish status (fishless versus fish bearing) in modulating historical food web response to dual forcing by bottom-up (resource regulation) and top-down (planktivory from stocked fish) processes in lakes of the southern interior of British Columbia. Specifically, we analyzed proxies of production in sediment cores over the past 100 years across multiple trophic levels including: (1) diatoms and diatom-inferred nutrient concentrations, (2) biochemical markers of algal biomass, and (3) sub-fossil Cladocera. Our study lakes have well-documented histories with respect to fisheries and watershed management. Two of our four study lakes were naturally fishless prior to stocking, whereas the other two lakes contained natural populations of native fish before being stocked. We hypothesized that if fish were important in structuring food webs in our lakes, we would document declines in the abundance of large daphnids and increases in biochemical markers of algal biomass (chlorophylls and carotenoids) post-stocking in our naturally fishless lakes, whereas little change would occur post-stocking in our naturally fish-bearing lakes.

STUDY AREA
The four study lakes are situated on the Thompson Plateau of south-central interior British Columbia and are within 40 km of the city of Kamloops along a north-south transect ( Figure 1). All lakes are situated within the dry and cool subzone of the Interior Douglas Fir (IDF) biogeoclimatic zone, between 940 and 1380 m elevation (Table 1). Mean annual air temperature (MAAT) and precipitation (MAP) recorded in Kamloops (elevation = 345 m) are 9. 4°C and 316 mm, respectively (1981-2010 average). Mean winter and summer temperatures are -1.6 and 20. 3°C, respectively (1981-2010 average) (Adjusted and Homogenized Canadian Climate Archive, http://ec.gc.ca/dcchaahccd). Although spatially variable, the MAP for the area near Roche Lake (Pratt Road, Figure 1) is 397 mm based on measurements from 1987 to 2006. Temperature recordings over this period suggest that mean summer and winter temperatures are 2.5 and 1.5°C cooler than in Kamloops, respectively, although this comparison is based only on intermittent recordings. Vegetation in our study area consists of stands of interior Douglas fir (Pseudotsuga menziesii), lodgepole pine (Pinus contorta), and spruce (Picea spp.), with some Trembling aspen (Populus tremuloides) and western redcedar (Thuja plicata). Geology of the Thompson Plateau is characterized by mid-Eocene and Miocene basalts that is mantled with glacial drift, with Precambrian granitic intrusions and batholiths occurring throughout the region, as well as numerous steepsided valleys, eskers, kames, and meltwater channels. Soils are predominantly gray luvisols and dystric brunisols (Valentine and others 1978). Figure 1. Location of the four study lakes in the southern interior region of British Columbia. The inset depicts the study region relative to the province of British Columbia. The shape of each lake basin is also shown, with associated coring location.  1957 1920-1930 1981-1990

STUDY LAKES AND MANAGEMENT HISTORY
Lakes were selected to have known histories related to fisheries management to be able to assess food web changes induced by stocking of salmonids. All study lakes are small (surface area £ 200 ha), moderately deep (10 m < Z max < 30 m), and meso-eutrophic and have moderate-to-high concentrations of dissolved organic carbon (Table 1). All lakes have been stocked quasi-annually predominantly with rainbow trout (Oncorhynchus mykiss) since the early to mid-1900s, with mean annual stocking densities of 200-300 fish ha -1 ( Table 1). Two of the lakes, Roche and Heffley, contained native fish populations prior to stocking and were treated with toxaphene a few years after stocking to eliminate native non-game fish species. Native non-game fish included: the northern pikeminnow (Ptychocheilus oregonensis), the peamouth chub (Mylocheilus caurinus), the redside shiner (Richardsonius balteatus), and the largescale sucker (Catostomus macrocheilus) (British Columbia Fisheries Inventory Data Queries, http://a100.gov.bc.ca/pub/fidq/). Based on fish surveys completed prior to stocking (for Peter Hope) and absence of chemical treatments to remove native non-game fish, both Peter Hope and Blackwell are believed to be naturally fishless prior to stocking (A. Klassen, pers. comm). Currently, all study lakes are monoculture systems of rainbow trout, with the exception of Peter Hope Lake where Kokanee (Oncorhynchus nerka) were introduced in 2017. Additional management activities over the past 100 years include the construction of outlet dams in the early 1900s on all lakes except Blackwell, and forest harvesting that has occurred within the watersheds at a range of intensities since 1970 (Table 1). Heffley Lake has 100 year-round residential structures around the shoreline, whereas Roche and Peter Hope lakes have 22 and 39 seasonal structures, respectively. Blackwell Lake has no structures in its watershed.

Sediment Core Collection
Sediment cores were collected in February 2016 from near the deepest portion of each lake using a Glew gravity corer (Glew and others 2001) with an internal diameter of 7.6 cm. All cores were sectioned at 0.25-cm intervals, and samples were placed into individual sterile bags. Sediments were transported to Queen's University (Kingston, Ontario) where they were stored in a dark cold room at 4°C until processed for analyses.

Chronologies
Chronologies for all sediment cores are based on measured activities of 210 Pb and 137 Cs via gamma spectroscopy that was carried out at the Paleoecological Environmental Assessment and Research Laboratory (PEARL), Queen's University, following procedures similar to Schelske and others (1994). Briefly, sediment samples were freeze-dried and packed into plastic tubes to 2 cm height, sealed with 2-Ton epoxy and left for at least 2 weeks to allow 214 Bi to equilibrate with 226 Ra. Radioactive decay was measured for 80,000 s using an EG&G Ortec germanium crystal well detector. Activities of total 210 Pb, 137 Cs and supported 210 Pb (via 214 Pb) were determined for all samples based on calibration to internal standards IAEA 312 and 385.
Two types of models were used to estimate dates for the intervals analyzed for radioisotope activities. For the cores from Peter Hope and Heffley lakes, a Constant Rate of Supply (CRS) model was used, determined from unsupported 210 Pb activities using the ScienTissiME program (Mike Scheer, unpublished). For the cores from Blackwell and Roche, a composite age-depth model similar to that presented in Appleby (2001) was used. These latter models use independent chronostratigraphic markers to estimate changes in 210 Pb supply rate when there are suspected discrepancies with a CRS model. We only had one reference point (that is, the 1963 CE 137 Cs peak associated with nuclear bomb testing), and thus, the CRS model was numerically constrained to this independent reference point following procedures of von Gunten and others (2009). This model was deemed appropriate for Blackwell due to a discrepancy of more than 10 years between the CRS estimated age and the well-resolved 1963 (CE) 137 Cs peak. Although the pattern of exponential decay of 210 Pb was not as strong in the Roche Lake core, an additional core with visible laminations collected and dated in 2018 showed a similar pattern in 210 Pb and 137 Cs activities. As such, we also used the constrained CRS model for Roche Lake. For all lakes, undated intervals were interpolated using a monotonic cubic spline (Simpson 2018a). All dates were expressed as years Common Era (CE).

Isotopes and Elemental Ratios
Stable isotopes of nitrogen (d 15 N) and carbon (d 13 C) were analyzed using a Thermoquest (Finnigan MAT) Delta Plus XL isotope ratio mass spectrometer (IRMS) coupled with a Carlo Erba NC2500 elemental analyzer following procedures Cladoceran Response to Fish Stocking and Nutrients of Savage and others (2010). Analysis was performed on freeze-dried sediments, and values were calibrated against the international standards of atmospheric N 2 and Vienna Pee Dee Belemnite (VPBD) for d 15 N and d 13 C, respectively. Isotopic data are presented as per mille (&). Sedimentary C/N is presented as mass ratios.

Chaoborus mandibles
Vertebrate planktivory often exerts a strong control over the abundance and composition of Chaoborus in lakes. In particular, C. americanus is considered to exist almost exclusively in fishless lakes due to their large size and non-migratory behavior (for example, Wissel and others 2003), whereas many other Chaoborus species are capable of vertical migration to avoid visually orienting predators. Therefore, we used the composition of sub-fossil Chaoborus mandibles prior to stocking in all lakes as an independent proxy of pre-stocking status of fish presence or absence. Chaoborus mandibles were isolated from wet sediment following protocols outlined in Walker (2001) and identified using the taxonomic key of Uutala (1990). If no mandibles were recovered after processing 2 g of wet sediment, analysis was discontinued.

Diatoms
For each sediment sample, a 50:50 molar solution of nitric (HNO 3 ) and sulfuric (H 2 SO 4 ) acid was added to 0.2 g of wet sediment in a 20-ml glass vial. Samples were then heated in a hot water bath at 75-80°C to help digest organic matter. After settling for 24 h, samples were aspirated to 5 ml and rinsed with double-deionized water. This procedure was repeated 7-8 times, until the samples reached the same pH as the deionized water. Samples were then aspirated to 5 ml, and a known volume of microsphere solution was added to each sample to enable the calculation of diatom concentrations. Four successive dilutions of each sample were then pipetted onto glass coverslips, allowed to air dry overnight, and then mounted onto glass slides using Naphrax Ò . At least 400 diatom valves were enumerated per sample, using a Leica DMRB microscope fitted with a 100 9 Fluotar objective using differential interference contrast at 1000 x magnification. Diatoms were identified to species level or lower, primarily using the taxonomic collections of Krammer and Lange-Bertalot (1986, 1988, 1991a, and Cumming and others (1995).

Pigments
Photosynthetic pigments (chlorophylls and carotenoids) were quantified in all cores to quantify historical changes in the abundance of algae and photosynthetic bacteria. Pigment analysis was performed using reversed-phase high-performance liquid chromatography (RP-HPLC) at the Institute of Environmental Change and Society (IECS), University of Regina, following procedures of Leavitt and Hodgson (2001). Briefly, pigments were extracted by soaking freeze-dried sediments in an 80:15:5 mixture (by volume) of acetone/ methanol/water for 24 h in low-light conditions. Pigment extracts were then dried under a stream of inert N 2 gas. An internal reference of Sudan II (3.2 mg l -1 ) was injected to each sample to resuspend pigment residues, which were then isolated from solution using an Aligent HPLC system calibrated with authentic standards from the Danish Hydraulic Institute (DHI) and from local pigment isolates (Leavitt and Hodgson 2001). Pigments were tentatively identified by chromatographic position and absorbance spectra, and concentrations were expressed as nmoles of pigment per gram organic matter (nmols g -1 organic matter).

Cladocera
Protocols outlined in Korhola and Rautio (2001) were followed to prepare samples for analysis of sub-fossil Cladocera. Briefly, 1 g of wet sediment was deflocculated in 150 ml of 10% KOH solution at 70°C for 30 min. Samples were then poured through a 34-lm sieve and rinsed with doubledeionized water. Material retained on the sieve was backwashed into a 12-ml vial, and a few drops of ethanol (preservative) and a safranin-glycerol mixture (dye) were added to each sample. Fifty microliter aliquots of each sample were pipetted onto glass slides and allowed to dry, and this process was repeated as necessary to achieve a sufficient amount of cladoceran remains for enumeration (minimum of 70 individuals; Kurek and others 2010). Glass coverslips were then mounted on the slides using glycerin jelly.
Remains were enumerated at either 200 x or 400 x magnification under bright-field illumination, using the taxonomic keys of Korosi and Smol (2012a, b) and Szeroczyñ ska and Sarmaja- Korjonen (2007). When possible, we also measured the length of the Bosmina carapace and attenule and the length of the post-abdominal claw from Daphnia pulex and Daphnia longispina, as changes in the size of these remains can provide additional insight into changes in predation pressure (Kitchell and Kitchell 1980).

Diatom-Inferred Total Phosphorus (DI-TP)
Total phosphorus (TP) was inferred through time in all cores from fossil diatoms (DI-TP) by application of a training dataset which related diatom species composition to measured TP values in 251 lakes in B.C. (Cumming and others 2015). Lakes in the calibration set used herein were located across several biogeoclimatic zones throughout B.C. yet were all freshwater (salinity < 500 mg l -1 ), and spanned a TP gradient from 2 to 227 lg l -1 . The model used to quantify changes in TP over time is based on TP optima estimates for individual taxa determined via weighted averaging regression on relative abundance data (Cumming and others 2015). The bootstrapped coefficient of determination (r 2 ) of the model is 0.53, with an RMSEP of 0.32. To be included in TP estimates, taxa had to reach at least 2% relative abundance. Additional information on the lake types included in the calibration set, model development, and model performance is described in detail in Cumming and others (2015).

Principal Curves and Generalized Additive Models
Principal curves (PrCs) were used to summarize compositional change through time for the diatoms in all cores, and the Cladocera in Roche and Heffley lakes. However, given that relatively little cladoceran compositional change occurred in Peter Hope and Blackwell lakes, we presented concentration data for these two systems. PrC is an ordinationbased method that fits a smooth curve through multivariate data, and is considered a form of nonlinear principal components analysis (De'ath 1999). PrCs were selected as they are considered to capture more variation in the multivariate data when compared to a PCA or CA, particularly when a single dominant gradient of compositional change is present (Simpson and Birks 2012). PCA axis-1 scores were used as starting points for curve fitting. Procedures outlined in Simpson and Birks (2012) were followed, allowing the complexity of smooth splines to vary, and using generalized cross-validation (GCV) to determine the spline degrees of freedom.
Given that paleoecological data are typically nonlinear and irregularly spaced in time (Simpson 2018a), we used generalized additive models (GAMs) to estimate trends in the sub-fossil cladoceran data and diatom data. This approach enabled us to assess whether historical changes in zooplankton coincided with initial fish stocking (topdown forcing), or with changes in diatom composition and associated TP estimates (bottom-up forcing). Diatom and cladoceran PrC scores were used as response variables, with the exception of Peter Hope and Blackwell, where cladoceran concentrations were modeled instead. In all cases, the temporal series of PrC scores (or concentration data) were modeled separately with 210 Pb age as the only covariate, using thin-plate regression splines as the basis spline, and restricted maximum likelihood (REML) to determine spline degrees of freedom. When cladoceran concentrations were modeled (as in Peter Hope and Blackwell lakes), a gamma error distribution was specified with a log link function. In all cases, initial models were fit with a continuous-time first-order autoregressive [CAR(1)] correlation structure (Simpson 2018a). However, in a number of instances this resulted in an unsatisfactory fit; in these cases, we followed the recommendations of Simpson (2018a) by increasing the size of the basis dimension and specifying observational weights for each sample to account for heteroscedasticity in the data. In the case where a model with observational weights was indistinguishable from one without weights, we selected the simplest model (that is, no observational weights). Periods of significant change were subsequently estimated by determining where the rate of change along the estimated trend was nonzero, as described in Simpson (2018a) and implemented in Bennion and others (2015). All statistical analyses were carried out in R using packages analogue (Simpson 2018b) and mgcv (Wood 2019).

Core Chronologies
With the exception of Roche Lake, radioisotopic analysis of sediment cores generally showed an exponential decline of unsupported 210 Pb with core depth, with an r 2 for the first-order exponential decay ranging from 0.90 to 0.98 ( Figure S1). Similarly, the r 2 of a linear model fit to the relationship between log 10 unsupported 210 Pb activity and cumulative dry mass ranged from 0.77 to 0.95. All lakes except Roche had well-resolved peaks in 137 Cs activity, serving as a strong independent marker for age. In Blackwell Lake, the 137 Cs-constrained CRS age-depth model was used to adjust the chronology for a discrepancy of more than 10 years between the 210 Pb-estimated age at 10 cm depth (ca. 1963 CE) and the 137 Cs peak ( Figure S1). Background (supported) 210 Pb was reached at core depths between 12 and 14 cm for Roche, Heffley, and Blackwell lakes, and 9 cm for Peter Hope Lake. However, the chronology for Roche Lake should be interpreted with caution, given the relatively poor fit of the exponential decay pattern (r 2 = 0.69) likely due to a large increase in sedimentation rate, and because we had to extrapolate the age-depth model to estimate the ages of pre-fish stocking intervals. Based on the dating models, sub-decadal temporal resolution was achieved for all proxies, on all cores, since 1900, with analysis undertaken at intervals ranging from 0.25 to 1 cm. Blackwell Lake was the only exception, with the core beginning at 1935. However, we note that this was sufficient to be able to characterize food web structure prior to fish stocking, which began in 1962 in Blackwell Lake.

Isotopes and Elemental Ratios
Patterns in C/N were broadly consistent among Roche, Heffley, and Blackwell lakes with progressive declines over the period of record. However, Roche and Heffley showed short-lived peaks in C/ N, occurring at ca. 1930-1940 CE, respectively, both at values of 13-14 ( Figure S2). Peter Hope also showed a broad period of elevated C/N, reaching a maximum of 20 by 1930 CE, after which it declined and remained relatively stable at a value of 13. d 13 C followed similar patterns as those exhibited by C/N in all cores, with correlations between the two time series varying between 29 and 93% of variance explained when performed on raw values and between 14 and 79% when performed on detrended (first differences) time series, suggesting terrestrial inputs have an important control on d 13 C in our cores. Carbon isotope values were similar (ca. -28 to -34&) in Heffley and Blackwell lakes, whereas more enriched values were recorded in the cores from both Roche and Peter Hope lakes, particularly during the mid-core peaks (ca. -20&) near 1930 ( Figure S2). In Roche and Peter Hope lakes, d 15 N followed a similar pattern as C/N with peaks at 5& ( Figure S2). Peter Hope Lake showed an additional period of elevated d 15 N ( 5&) post-1980, whereas d 15 N remained relatively stable at 4& for the remainder of the record in Roche Lake. Blackwell Lake showed a progressive decline in d 15 N similar to that of the C/ N, changing from 3& in the early 1900s, to less than 2& by 1990, whereas a progressive enrichment of d 15 N by 2& occurred in Heffley Lake between 1900 and 1990 CE ( Figure S2).

Ecological Change
Identification of Chaoborus mandibles is consistent with Peter Hope being fishless prior to stocking, with elevated abundance of C. americanus near the base of the core (Table S1). Although 14% of the pre-stocking chaoboriid assemblage in Blackwell Lake was composed of C. americanus, C. trivittatus was the predominant taxon. This latter species is also considered sensitive to fish predation due to its weak vertical migration and large size (Lamontagne and Schindler 1994). In contrast, C. flavicans and C. sayomyia were present in Roche Lake, which are common in food webs with zooplanktivorous fishes (Uutala 1990;Lamontagne and Schindler 1994). No mandibles were recovered from pre-stocking sediments in Heffley Lake (Table S1). Moreover, given that native non-game fish were caught in Roche and Heffley lakes around the time of stocking, we contend that these basins supported native fish populations. These lines of evidence support our study design based on native fish status and allow us to interpret subsequent trophic changes within this context.

Historically Fish-Bearing Lakes (Roche and Heffley)
A total of 143 and 104 diatom species were identified in the cores from Roche and Heffley lakes, respectively, although only 20% and 25% of these taxa reached relative abundances above 2% in at least one sample in Roche Lake and Heffley Lake, respectively.  S4). Cores from both lakes showed an abrupt switch from an assemblage dominated by oligo-to mesotrophic diatom taxa (for example, C. bodanica v. lemanica, L. michiganiana, F. crotonensis, A. formosa) to an assemblage that consisted of more than 80% of the more eutrophic Stephanodiscus taxa in the mid-1900s (Figures S3 and S4). This change was reflected in the DI-TP, increasing from a mean of 16-34 lg l -1 after ca. 1920 in the core from Roche Lake, and a mean of 12-34 lg l -1 after ca. 1940 in the core from Heffley Lake (Figure 2A, B). In both lakes, this shift to a more nutrient-rich state was marked by increases in total diatom concentration ( Figures S3 and S4). Roche Lake also showed an additional increase in DI-TP after ca.
2000, which was driven by increases in both S. hantzschii and S. hantzschii tenuis. GAMs fitted to diatom PrC scores from Roche Lake showed three periods of significant compositional change, between ca. 1900-1925, 1940-1970, and after 2000, which corresponded to the DI-TP (Figure 2A). Similarly, the core from Heffley Lake showed sig- Figure 2. Summarized paleolimnological data from multiple trophic levels for the two naturally fish-bearing lakes A Roche Lake and B Heffley Lake. The solid line indicates the approximate date of outlet dam construction on Roche Lake. The first dashed line (fish symbol) indicates date of initial stocking, while the second dashed line (triangle) indicates date of first toxaphene treatment. Bold regions on generalized additive models (GAMs) fitted to the cladoceran and diatom data indicate periods of significant change. Algal affinities for each pigment are indicated in italic and are based on Leavitt and Hodgson (2001). All pigments are in units of nmols g -1 organic matter. Although Heffley Lake does have an outlet dam, we could only constrain its date of construction between 1925 and 1946. nificant compositional change beginning in ca. 1930, a period which also coincided with the increase in DI-TP ( Figure 2B).
Concentrations of pigments representing the biomass of important phytoplankton groups including diatoms (diatoxanthin), cryptophytes (alloxanthin) and total cyanobacteria (echinenone) showed a strong response to elevated DI-TP in Heffley Lake, with most biomarkers exhibiting a two to threefold increase in concentrations in the mid-1900s ( Figure 2B). Similar patterns were recorded for ''bloom-forming'' phytoplankton derived from chlorophytes and cyanobacteria (luteinzeaxanthin) and total phototrophic abundance (as b-carotene). However, unlike the trend in DI-TP, most of these pigments declined to concentrations similar to those prior to eutrophication, with the exception of stable canthaxanthin from colonial Nostocales cyanobacteria and labile fucoxanthin from siliceous algae and some dinoflagellates (Figures 2B, S4). Okenone from obligately anaerobic purple sulfur photosynthetic bacteria was recorded throughout the core, suggesting at least seasonal anoxic habitats were present over the past ca. 100 years.
Pigments in Roche Lake sediments showed subtle changes, with a period of elevated concentrations of biomarkers from diatoms (diatoxanthin), cryptophytes (alloxanthin), chlorophytes + cyanobacteria (lutein-zeaxanthin), and total cyanobacteria (echinenone) in the early 1900s that coincided with the first period of significant compositional change in the diatoms, and a small peak in DI-TP ( Figure 2B). However, these pigments recovered by ca. 1950, before exhibiting an abrupt increase after ca. 2000, coeval with marked increases in DI-TP. As in Heffley Lake, canthaxanthin from Nostocales cyanobacteria did not record the initial transient peak, and instead, along with labile fucoxanthin, increased markedly only after ca. 2000. General pigments from colonial cyanobacteria (myxoxanthophyll), as well as compounds characteristic of potentially N 2 -fixing forms (aphanizophyll) were also present in Roche Lake, and increased markedly after ca. 2000 ( Figure S3). As in Heffley Lake, okenone from purple sulfur bacteria were present throughout core ( Figure S3). In cores from both systems, the ratio of labile chlorophyll a to its stable degradation product pheophytin a, an index of preservation (Leavitt and Hodgson 2001), showed little change through time (Figures S3 and  S4).
In both lakes, pelagic cladocerans including Bosmina spp. and members of the Daphnia pulex and Daphnia longispina complexes were most dominant, with maximum abundances ranging from 40 to 90% of the fossil assemblage ( Figures S3, S4). Littoral taxa including Chydorus brevilabris and Alona spp. were also present throughout the cores, albeit at low abundances (< 15%). Alona consisted predominantly of A. circumfibriata/A. guttata (which are difficult to differentiate based on headshields alone (see Korosi and Smol 2012b), A. quadrangularis, and A. affinis. Both lakes showed abrupt shifts in Cladocera species composition. Initially, assemblages were dominated by Bosmina spp. (> 50%) and D. longispina complex prior to changing to an assemblage dominated by large-bodied D. pulex in both Roche (ca. 1920;Figures 2A, S3) and Heffley lakes (ca. 1935;Figures 2B, S4). GAMs fit to PrC scores showed that these changes in composition were significant (Figure 2A, B). In both lakes, cladoceran changes coincided closely with the increases in diatom-inferred TP, as well as the date of initial stocking in Heffley Lake. Analysis of concentration data suggests both declines in Bosmina spp. and increases in the daphnids (Figures S3, S4).

Historically Fishless Lakes (Peter Hope and Blackwell)
A total of 96 and 147 diatom taxa were identified in the cores from Peter Hope and Blackwell lakes, respectively, but only 30% of these taxa reached relative abundances above 2% in Peter Hope Lake, and 20% in Blackwell Lake. However, these taxa represented at least 80% of the assemblage in both basins. Important planktonic taxa throughout these cores were similar to those in Roche and Heffley, although Blackwell Lake had some additional taxa including Aulacoseira (subarctica (Otto Mü ller) E.Y. Haworth, and valida (Grunow) Krammer), Discostella pseudostelligera (Hustedt) Houk &Klee, Fragilaria capucina v. gracilis (Østrup) Hustedt, and Tabellaria flocculosa str 3p (Roth) Kü tzing. Both lakes showed periods of increased DI-TP, with Peter Hope Lake increasing from 20 to 30 lg l -1 at ca. 1930 and Blackwell Lake increasing from 20 lg l -1 in the early 1900s to a maximum of 35 lg l -1 in the early 1980s ( Figure 3A, B). These changes were primarily driven by an increase in the proportion of meso-eutrophic Stephanodiscus taxa, and A. formosa to some extent in the core from Blackwell Lake ( Figures S5, S6). Unlike the trends in the cores from both Roche and Heffley lakes, DI-TP in both of these systems recovered to inferred TP concentrations similar to those in the early 1900s in both Peter Hope (by ca. 1970) and Blackwell lakes (by ca. 2000). Notably, the period of elevated DI-TP in Peter Hope Lake coincided with an increase in diatom concentration of two orders of magnitude, relative to the concentrations between 1900 and 1920 ( Figure S5). Blackwell Lake showed some additional changes in species composition, with a post-1990 increase in a number of benthic taxa including Staurosirella pinnata (Ehrenberg), Pseu-dostaurosira brevistriata (Grunow), and Navicula spp. ( Figure S6). GAMs fitted to PrC scores showed significant compositional changes that followed the trend in DI-TP in the Peter Hope core ( Figure 3A). Significant compositional change occurred in the Blackwell Lake core that followed the initial trend Figure 3. Summarized paleolimnological data from multiple trophic levels for the two naturally fishless lakes A Peter Hope Lake and B Blackwell Lake. The solid line indicates the approximate date of outlet dam construction on Peter Hope Lake. The dashed line (fish symbol) indicates date of initial stocking. Bold regions on generalized additive models (GAMs) fitted to the cladoceran and diatom data indicate periods of significant change. Units for D. pulex concentrations are 9 10 3 g -1 dry wt. for Peter Hope Lake and 9 10 2 g -1 dry wt. for Blackwell Lake. D. pulex PCL represents the average length of the post-abdominal claw. Algal affinities for each pigment are indicated in italic and are based on Leavitt and Hodgson (2001). All pigments are in units of nmols g -1 organic matter.
in increasing DI-TP in the mid-1900s. Trends in DI-TP and GAMs fitted to PrC scores diverged at 1980, as a unique assemblage of more mesotrophic taxa (F. crotonensis, T. flocculosa str 3p), and a number of benthic taxa (P. brevistriata, S. pinnata, and Navicula spp.) increased in abundance (Figures 3B, S6).
Analysis of fossil pigments in the core from Peter Hope Lake reflected a change to more eutrophic conditions (inferred from the diatoms), with increases in diatoms (diatoxanthin), cryptophytes (alloxanthin), chlorophytes + cyanobacteria (lutein-zeaxanthin), and Nostocales colonial cyanobacteria (canthxanthin) at ca. 1930 (Figure 3A). A number of these pigments showed subsequent declines that coincided with reduced DI-TP, stabilizing at lower concentrations by ca. 1970, although not as low as pre-1930 values. The core from Peter Hope Lake also included some additional diatom-related pigments, including diadinoxanthin and fucoxanthin which both followed patterns similar to diatoxanthin ( Figure S5). Okenone was also present in Peter Hope Lake after 1950.
Concentrations of most pigments in the core from Blackwell Lake showed a pattern different from the other study lakes, with declines in compounds from diatoms (diatoxanthin), cryptophytes (alloxanthin), and chlorophytes + cyanobacteria (lutein-zeaxanthin) that occurred even when DI-TP was high ( Figure 3B). However, subtle increases in markers for colonial cyanobacteria (canthaxanthin) and total cyanobacteria (echinenone) did occur and coincided with elevated DI-TP between ca. 1960 and 2000 ( Figure 3B). Lastly, in cores from both systems, the ratio of chlorophyll a to pheophytin a showed little change through time ( Figures S5, S6), suggesting a consistent preservation environment.
Cores from Peter Hope and Blackwell lakes showed little change in cladoceran species composition through time and were dominated by large pelagic Cladocera of the D. pulex complex (> 50%), and D. longispina to a lesser extent (Figures S5, S6). However, significant increases in the concentrations of D. pulex occurred in both cores, which coincided closely with diatom assemblage change and the timing of increased DI-TP ( Figure 3A, B). These changes occurred shortly after initial stocking in both lakes. We also note that in contrast to the cores from Roche and Heffley lakes, fossil assemblages in these two naturally fishless lakes were dominated by large daphnids, with little to no Bosmina, both before and after stocking.

DISCUSSION
Traditional food web theory predicts that in the absence of strong vertebrate planktivory, large zooplankton will outcompete small zooplankton because they are more effective grazers, whereas small zooplankton will dominate when vertebrate planktivores are present due to size-selective predation on large taxa (Carpenter and others 1985). The native cladoceran assemblages prior to fish stocking in our study lakes were consistent with our a priori-defined groups with respect to food web status (that is, naturally fishless or naturally fish bearing) with mainly large-bodied Daphnia spp. and C. americanus and C. trivittatus in fishless Peter Hope and Blackwell lakes, and abundant smallbodied Bosmina spp. in naturally fish-bearing Roche and Heffley lakes. Therefore, within the context of our study design, declines in the large daphnids would be expected to occur post-stocking in Peter Hope and Blackwell, while little change would be expected in Roche and Heffley given that cladoceran assemblages were pre-shaped by native fish. Contrary to our prediction, the concentrations of large daphnids (primarily D. pulex) increased following stocking in three of four lakes and remained high following stocking in one lake (Roche Lake). Because the Pennask strain of rainbow trout that are most commonly stocked in our lakes are known to be predominantly pelagic or mid-water foragers (FFSBC 2004) and often feed on Daphnia (Hirsch and Negus 2000; Watson and others 2019), it is surprising to find no evidence of a change in predation regime in our naturally fishless lakes, especially given that these lakes are stocked at densities that have induced shifts in zooplankton size structure in other systems (Leavitt and others 1994;Donald and others 2001). The absence of changes in zooplankton size structure following stocking was also supported by the size measurements of the D. pulex post-abdominal claw length (PCL) in our two fishless lakes, which showed little change ( Figure 3A, B). Instead, all of our study lakes experienced sustained or episodic periods of eutrophication over the past 100 years. Given the coinciding trends between the onset of increased DI-TP, algal abundance (from pigments) and D. pulex abundance in our lakes, we infer that direct or indirect effects of eutrophication led to the observed changes in zooplankton abundance and community composition. Specifically, we suggest that the intensification of deepwater anoxia associated with eutrophication may have provided additional refugia for D. pulex from visually orienting predators, forcing rainbow trout to utilize alternative food sources. Moreover, because Daphnia may face reproductive costs associated with inhabiting cold, anoxic refugia, we also infer that increased primary production favored increased secondary production of Daphnia.

Post-stocking Daphnia Success and Changes in Predator Refugia
The post-stocking increases in Daphnia relative abundance and concentration in our study lakes indicate that their populations have not been diminished by predation from rainbow trout, and we infer that this arises due to an interaction between changes in the availability or extent of refugia from zooplanktivores, and the basic trophic ecology of the introduced rainbow trout. Deepwater anoxia, or enhanced severity of anoxia, often accompanies eutrophication. Biweekly monitoring data confirm that severe seasonal hypolimnetic oxygen depletion occurs in both Roche and Peter Hope lakes (no data for other sites), with dissolved oxygen concentrations of 1-2 mg l -1 below 8 m in the water column in Roche Lake and below 16 m in Peter Hope Lake ( Figure S7). Recent work suggests that some large-bodied daphnids can utilize the anoxic zone as a refugium from predators. For example, Larsson and Lampert (2011) showed that mortality of Daphnia pulicaria (considered part of the D. pulex complex; Hebert 1995) from predation was reduced by over 50% when provided a hypoxic hypolimnion in a large mesocosm experiment. This response was attributed to that fact that Daphnia can better withstand oxygen depletion by synthesizing hemoglobin (Weider and Lampert 1985). Although rainbow trout can tolerate concentrations of dissolved oxygen considered marginal for other salmonids (for example, lake trout), they still prefer well-oxygenated waters, with physiological and behavioral response to oxygen stress occurring at 5 mg l -1 (Kerr and Lasenby 2000). This observation suggests that in our lakes, Daphnia may be utilizing a larger hypoxic zone to avoid predation. Hembre and Megard (2003) provided field evidence for this phenomenon in a small Minnesotan lake containing rainbow trout, where large D. pulicaria exploited a deep hypoxic (O 2 3-5 mg l -1 ) zone from June to October. Indeed, comparison of stomach contents of rainbow trout in nutrient-rich lakes in our study region with and without deepwater refuges support our inference that they are important in limiting fish predation on Daphnia (Tattersfield 2006). Although the presence of okenone in cores from Roche and Heffley lakes suggests that these basins were predisposed to provide hypoxic deepwater refugia for zooplankton prior to nutrient enrichment, the volumetric extent has presumably intensified over the past 100 years such that it serves as a more effective refuge with the introduction of a less hypoxia-tolerant salmonid.
Shifts in invertebrate predation may have also played a role in reshaping cladoceran assemblages, particularly in our naturally fish-bearing lakes. Previous studies have shown that rainbow trout may spend a high proportion of time in littoral and offshore epilimnetic areas (for example, Warner and Quinn 1995), which may be intensified in our lakes due to enhanced seasonal anoxia. Consequent reductions in vertebrate planktivory within pelagic habitats may have favored development of invertebrate predators. In particular, migratory Chaoborus (that is, not C. americanus) tend to feed on small zooplankton, particularly Bosmina (for example, Sutor and others 2001) and may be high enough to limit Bosmina production in some cases (Yan and others 1991). Such a food web change may have contributed to the decline in density of Bosmina in Roche and Heffley lakes. Further research, including a more comprehensive analysis of chaoborid fossils in additional cores, will be required to evaluate any potential cascading influence of limited vertebrate planktivory in pelagic areas on invertebrate predators and small herbivorous zooplankton.
Lateral migration of zooplankton may have also altered daphnid susceptibility to zooplanktivorous trout. Some studies have documented diel horizontal migration (DHM) of Daphnia into littoral areas during daytime to avoid visual predators (Stansfield and others 1997). In theory, eutrophication and potential expansion of littoral vegetation could provide Daphnia with additional predator refuge. However, DHM has largely only been detected in shallow, polymictic systems where vertical migration may not be advantageous (Burks and others 2002). In contrast, our lakes are relatively deep ( 20-30 m with the exception of Blackwell), and available monitoring data for Roche and Peter Hope lakes show marked seasonal thermal stratification ( Figure S7). Moreover, with the exception of Blackwell, concentrations of benthic diatoms exhibited either declines or little change when our study lakes became more nutrient rich ( Figures S3-S6), suggesting that benthic and littoral habitat has not expanded over time. Nevertheless, it has been suggested that if alternative benthic prey is highly abundant for vertebrate predators, littoral regions may still serve as an effective refuge for Daphnia (Burks and others 2002). Therefore, although there is little direct evidence for DHM in our lakes, the possible expansion of the importance of littoral areas for Daphnia requires further experimental work.

Daphnia Production and Resource Availability
Changes in food quantity and quality are important factors in the production of zooplankton. Although deepwater refugia is likely important for Daphnia survival, there are also reproductive costs associated with inhabiting cooler anoxic or hypoxic refugia (Hanazato and others 1989;Larsson and Lampert 2011). Because not only the relative abundance of Daphnia increased, but also total zooplankton population (as fossil concentration), it suggests that bottom-up factors may have been important. Indeed, analysis of photosynthetic pigments showed that the abundance of many algal groups increased as DI-TP increased, although this pattern was not as pronounced in Blackwell Lake (Figures 2, 3, S3-S6). Unicellular and flagellate groups of algae including diatoms, cryptophytes, dinoflagellates, and chlorophytes increased, which are viable food sources for zooplankton, including D. pulex (that is, food particle size £ 40 lm; Burns 1968). Daphnia are effective herbivores and, with the decline in pelagic Bosmina in Roche and Heffley lakes, would have less competition for algal resources (Vanni 1986). Such bottom-up regulation of Cladocera may have been particularly important in our two naturally fishless lakes, given that it is not possible that predation from fish was weakened following the introduction of rainbow trout.
Constant or increasing abundance of colonial cyanobacteria (as canthaxanthin, myxoxanthophyll, aphanizophyll) suggests that eutrophication did not greatly improve food quality. These cyanobacterial pigments showed prolonged or episodic increases associated with elevated inferred nutrients in our study lakes, particularly in Roche Lake after ca. 1990 (Figures 2, 3, S3-S6). In highly eutrophic lakes where cyanobacteria are abundant, Bosmina often dominates because of its ability to feed more selectively (DeMott 1982), whereas the less-selective Daphnia may suffer declines due to consumption of nutritionally inadequate colonial cyanobacteria, or mechanical interference of filtering by thoracic limbs (de Bernardi and Giussani 1990;DeMott and others 2001). Given that we did not observe declines in Daphnia when cyanobacteria were abundant, we infer that Daphnia maintained a competitive advantage in these lakes during productive periods, potentially augmented by deepwater refugia.

Toxaphene
Both Roche and Heffley lakes were treated with toxaphene 10 years following the initial introduction of rainbow trout to eliminate native nongame fish and reduce interspecific competition for resources (Table 1). Such treatments can impact invertebrate communities. For example, Miskimmin and  attributed declines in sedimentary Bosmina remains to toxaphene treatment in a small eutrophic lake in Alberta. Neeman (1962) provides one of the few monitoring records of zooplankton response to whole-lake toxaphene treatment, and in all four lakes studied, Bosmina remained stable approximately 1 week post-treatment, before showing declines in abundance of varying degrees up to 1-year post-treatment. This suggests that Bosmina reductions may have been a result of reduced vertebrate planktivory on large Cladocera, rather than a toxic effect per se. Moreover, the toxaphene concentration used to treat Heffley Lake did not exceed the recommended non-lethal limit for Bosmina (0.01 ppm). This, combined with the observation that neither Peter Hope nor Blackwell lakes were ever chemically treated, suggests that toxaphene toxicity cannot explain the common pattern of increased daphnid density observed in all study lakes (Figures 3, S3-S4).

Potential Drivers of Eutrophication
The introduction of game fish to small lakes may have important implications for nutrient cycling in small lakes. For example, mesocosm and wholelake studies have demonstrated that fish may increase limnetic P availability (Vanni and others 1997;Schindler and others 2001), a pattern that may be attributed to both nutrient excretion and remineralization of fish biomass. In our lakes, fish stocking coincided closely with increased DI-TP (Figures 2, 3), indicating it may have been an important source of nutrients. However, even under the assumption that the annual addition of fish biomass to each lake was completely remineralized, P associated with fish accounted for only a small fraction ( £ 2%) of the phosphorus needed to raise lakes from the minimum DI-TP pre-stocking to the maximum DI-TP observed post-stocking (Table S2), even assuming a water residence time 3-4 9 longer than those estimated for our lakes. Thus, although already a small portion of P, these values likely represent overestimates of the role that fish addi-tion played in the transition to more nutrient-rich systems.
The significant changes in primary and secondary producers in the early 1900s also coincided with the construction of outlet dams on Roche and Peter Hope lakes (Figures 2A, 3A), and potentially Heffley Lake ( Figure S3). Lake impoundment can raise water levels and inundate new terrestrial and littoral areas leading to a trophic surge for decades after impoundment (Thornton and others 1990). Consistent with this mechanism, all lakes except Blackwell (never dammed) showed increases in C/ N around the estimated time of damming, a pattern characteristic of enhanced terrestrial inputs (Meyers and Ishiwatari 1993). Similarly, the abrupt increase in d 15 N in Roche Lake ( Figure S2) may have been caused by leaching of N-rich soil following impoundment (for example, Das and others 2009). However, the lack of recovery in DI-TP in Roche and Heffley lakes is inconsistent with the theory of reservoir ontogeny (that is, trophic surge followed by recovery; Thornton and others 1990), possibly reflecting continued internal P loading due to seasonal anoxia augmented by recent changes in climate, or continued influx of allochthonous nutrients from shoreline development or logging (see site description). Indeed, previous paleolimnological assessments of eutrophication in interior B.C. lakes with human disturbance have shown similar trends in diatom species composition and associated changes in TP estimates over the last century (Reavie and others 2000). Lastly, while Blackwell Lake showed a period of elevated DI-TP throughout the mid-to late 1900s ( Figure 3B) but was never impounded, it did coincide closely with a period of high MAP from 1950 to the early 2000s ( Figure S8), which similarly may have increased the delivery of nutrients from the watershed. Given that Blackwell Lake has experienced the least amount of anthropogenic impact (Table 1), climate-mediated changes in the delivery of particles and solutes may be more pronounced for this system (Leavitt and others 2009).

CONCLUSIONS
This study explored the relative importance of historical nutrient and vertebrate planktivore controls on cladoceran assemblages in four productive fisheries lakes in interior British Columbia that have been stocked with rainbow trout since the early twentieth century. We documented increased abundance of large Daphnia post-stocking in all lakes, regardless of pre-stocking fish status (that is, fishless or fish bearing), a pattern that is opposite to expectations associated with changes in size-selective planktivory by fishes. Although vertebrate planktivory is often responsible for declines in large zooplankton following stocking, we suggest that an interaction between: (1) the development of deepwater anoxic refugia for Daphnia associated with eutrophication, (2) the introduction of an anoxia-intolerant planktivore (rainbow trout), and (3) increased resource availability to zooplankton, were all important drivers of post-stocking Daphnia success. Therefore, we suggest that both direct and indirect effects of bottom-up forces (primary production, anoxia, respectively) have enabled Daphnia to thrive, but these conditions were also a prerequisite for limiting vertebrate planktivory by trout. Given the apparent difference in the processes that altered food web structure in our lowerelevation lakes compared to high-elevation mountain lakes (for example, Schindler and others 2001), our study provides important insights on how lake-type impacts response to stocking, which may be important for managers who develop or maintain stocking programs.