Magnitude and regulation of zooplankton community production across boreal lakes

A major outstanding question in plankton ecology is whether the regulation of zooplankton production at the community level follows the same patterns that have been observed for individual populations. We used a novel biochemical approach to estimate in situ rates of crustacean zooplankton community production in 83 boreal lakes, with the objective of identifying the main drivers of zooplankton production at the community level across the boreal landscape. Our results show that the relationship of zooplankton community production to average community body size, total biomass, and temperature is comparable to what has been observed for individual populations. At the community level, however, there are additional drivers, including lake morphometry, and catchment properties, which further influence zooplankton production and which cannot be inferred from population level patterns.

The transfer of energy and biomass from primary producers to higher trophic levels is central to the functioning of aquatic ecosystems. Zooplankton are the major trophic intermediary in lakes, and understanding the dynamics of material and energy transfer through aquatic food webs necessarily requires quantification of biomass production rates of zooplankton communities, and an understanding of factors regulating these rates. Notably, whereas both the magnitude and drivers of phytoplankton and fish production at the community level in north-temperate lakes have been extensively studied (Schindler 1974;Vollenweider 1976;Karlsson et al. 2009;Benoît et al. 2016), the magnitude and regulation of zooplankton community production have not. Traditional approaches to quantify zooplankton production are based on extremely labor-intensive and detailed temporal monitoring of cohorts of individual populations (Downing and Rigler 1984;Dolbeth and Cusson 2012). Alternatively, approaches relying on biomass size spectra analyses (Sprules and Munawar 1986;Sprules and Stockwell 1995) can provide "community-level" information on productivity and transfer efficiency, but these analyses often require prey size measurements and assume fixed (often relative) transfer efficiencies.
Cohort-based approaches yield estimates of production that apply only to the target populations, yet it is community production which is at the core of many ecological and trophic questions (e.g., Garcia-Comas et al. 2016). For instance, quantitative evaluation of trophic interactions using productivity estimates for bacterial and primary producers (typically community-level measurements) require similar communitylevel zooplankton productivity measurements. Community production must be estimated individually on the basis of cohort analysis, rendering this population-based approach to community production difficult if not impossible to implement across many lakes, especially in remote regions where repeated sampling is a challenge. It is for this reason that other studies have not measured actual population production rates, but instead have used available empirical models to estimate community production (Johannsson et al. 2000;Stewart et al. 2010;Kelly et al. 2014;Mehner et al. 2015) based on published models linking population production to population biomass, body size, and water temperature (Downing 1984;Shuter and Ing 1997;Stockwell and Johannsson 1997).
Whereas population-based approaches such as summed population estimates or population-based empirical models may provide reasonable estimates of community production, neither approach lends itself to exploring the regulation of zooplankton production at the community level across lakes and along environmental gradients. On one hand, it is difficult, if not impossible, to infer the factors regulating zooplankton production at the community-level on the basis of regulation of the individual populations composing the community, because at the community-level there may be additional drivers, such as those determining the partition of total biomass into the different populations, and ecological interactions between constituent populations (Gliwicz 2009). On the other hand, it is unclear whether zooplankton production at the community-level follows scaling patterns relative to total biomass, mean body-size, and temperature resembling those identified for the constituent populations. Plante and Downing (1989) demonstrated consistency in both the body size and temperature dependencies of production across a wide range of invertebrate species, suggesting that a mixed community may indeed follow the general patterns of its component individual populations. This conclusion, however, is based on aggregating data across species across many orders of magnitude in body size, from rotifers to large benthic invertebrates; the pattern (Plante and Downing 1989) appears to break down within narrower ranges of body mass. In particular, there is evidence that for lake zooplankton, specific production (or production to biomass ratio; P/B) varies greatly across species (Banse and Mosher 1980), and also through time and space within a single species (Downing 1984). The extent to which population-level metabolic scaling can be extrapolated to the community level is a question that is relevant not only to our understanding of lake plankton, but to other communities as well, including phytoplankton, invertebrates, and fish, and goes to the core of ecology, for example, underpinning the metabolic theory of ecology (Brown et al. 2004).
Assessing the regulation and body-size scaling of zooplankton production at the community level therefore requires an independent quantification of total production rate that is not based on population-level estimates, and this across a variety of lakes and along sufficiently broad environmental gradients. Here, we present a large-scale study of the magnitude and regulation of zooplankton community production, using a novel approach to quantify community-level crustacean zooplankton biomass production rates that does not involve extrapolating population-based empirical models. Our method is based solely on measurement of the standing activity and rate of decay of the crustacean moulting enzyme, chitobiase , which yields an integrative, short time scale, measure of community production rates that can be obtained relatively easily across many lakes. We applied this approach across a large suite of boreal lakes that span wide geographic and environmental gradients to determine: (1) the range in daily total crustacean production rates across lake gradients; (2) the mean body-size and temperature dependency of community production rates in comparison with those reported for individual zooplankton populations; (3) the main environmental drivers of short time scale zooplankton community production, and the most robust predictive models. We use these elements to develop an integrative scheme for zooplankton community production regulation across boreal lakes, combining elements of zooplankton community structure and key environmental factors at the lake and catchment scales.

Study lakes and sampling
Daily total crustacean zooplankton production rates (ZP) were measured in 83 boreal lakes from three regions across Quebec, Canada (Supporting Information Fig. S1, Table 1). Twelve lakes in the Chibougamau region (49.10-50.858N, 72.98-74.528W), and 27 lakes in the Chicoutimi region (47.78-48.778N, 71.13-72.128W) were sampled from 03 June 2011 to 02 October 2011. Some Chibougamau and Chicoutimi lakes were sampled up to three times (five lakes for three times and seven lakes for two times), with the temporal separation between these being sufficiently long (38-74 d) to assume novel environmental conditions and communities composed of different species or population cohorts at very different stages of seasonal development, and thus independent daily production estimates. A further 43 lakes were each sampled once in the Schefferville region (54.16-55.478N, 65.21-68.848W) in August Daily biomass production rates were estimated from measurements of the standing native activity (CBA NAT ), and the measured turnover rate (T CBA ) of the crustacean molting enzyme, chitobiase (CBA), in the water column using samples collected at 2 m depth (following Sastri et al. 2013). We assumed that samples from 2-m depth captures both the portion of the community that carries out diurnal migration but which nevertheless leaves a moulting trace in the epilimnion, as well as the portion of the community that resides predominantly in the epilimnion, as observed elsewhere for similar lakes (Brosseau et al. 2012); these estimates may nevertheless represent underestimates, especially in terms of migrating Chaoborus, which inhabit bottom waters and are seldom present in surface water during daytime, when samples were taken. Daily ZP represents the production rate of biomass attributable to all somatic growth for the entire developing (actively molting) crustacean zooplankton community (methods in Supporting Information). Total zooplankton biomass across the 105-2000 lm equivalent spherical diameter size range was calculated from Laser Optical Plankton Counter measurements for each net sample (see Supporting Information). Mean individual body mass for the community was calculated as the product of abundance and the mid-point of mass within a size class summed across all size classes > 90 lm and divided by total abundance.
To explain cross-lake variation in ZP, we measured a suite of environmental variables (within-lake and catchment levels, Table 1). Mean epilimnetic temperature (8C) was estimated from a multi-parameter sonde (YSI, Yellow Springs Instruments) profiles. Samples for phytoplankton biomass (chlorophyll a [Chl a], lg L 21 ), total phosphorus (TP, lg L 21 ), and dissolved organic carbon (DOC, mg L 21 ) were collected from 0.5 m depth. In the lab, Chl a was extracted with 90% hot ethanol and measured spectrophotometrically before and after acidification (after Lorenzen 1967;Nush 1980); TP was measured using the molybdenum-blue method following persulfate digestion (Cattaneo and Prairie 1995); and DOC concentration was measured on an O.I. Analytical TIC/TOC using 0.45 lm filtered water after sodium persulfate digestion. Lake area was derived using ArcGIS V10 and the National Topographic Data Base and estimated catchment forest and wetland cover obtained from Geobase (2009).

Statistical analyses
To compare patterns in total daily crustacean production rates with those reported for individual populations, and in particular with the Plante and Downing (1989) model based on a compilation of species-specific production data, we first built a regression model containing similar variables: standing biomass, mean size, and temperature (Table 2, model A). We also explored potential environmental drivers of ZP, using "best subset regression" to identify the best models for a selected number of explanatory variables (Table 2, models B-G) based on the significance (p < 0.05) of both the overall regression, ranked according to akaike information criterion (AIC) values, and the explanatory variables. Finally, we merged the community structure, environmental, and catchment-level variables (Table 2, model H). This sequential approach allowed the emergence of the best combined predictive ZP model.
We further explored the relationships between the various drivers of ZP using structural equation modeling (SEM ;Bollen 1989;Shipley 2002;Grace 2006). This approach allowed us to: (1) assess whether environmental variables affected ZP directly via growth rates, or indirectly via community structure (biomass and body size); and (2) to better delineate the complex relationships between environmental factors and ZP (Fig. 1a). factors that directly or indirectly affect zooplankton production in lakes (Fig. 1a). To find the most parsimonious SEM that could explain our data, we used an iterative fitting procedure. At each step, we either added or removed the link that resulted in the largest decrease in the overall model AIC with no effect on the R 2 of ZP. The fit of each SEM to our data was assessed using the comparative fit index (CFI) and the p-value of the robust chi-square (v 2 ), where a significant SEM model is one having a p-value greater than 0.05, as the null hypothesis is that there are no significant differences between the SEM and the data. Variables with no direct or indirect path to biomass   Table S1). Variables are: Area, lake area; Biom, standing biomass; Chl a, chlorophyll a; DOC, dissolved organic carbon; %FOR, percent forest cover in catchment; Size, mean individual body mass; Temp, epilimnion temperature; TP, total phosphorus; %WET, percent wetland cover in catchment; ZP, zooplankton biomass production.
or ZP were removed from the final SEM. SEM analyses were carried out with the Lavaan package (Rosseel 2012) in R version 3.03 (R Core Team 2013).

Results
Lakes varied widely in size, trophic status, and temperature (Table 1). Zooplankton community biomass across all lakes ranged more than two orders of magnitude (3.3-1147 mg C m 23 ). Mean individual body mass across the crustacean communities ranged from 0.17 lg to 3.36 lg per individual, while chitobiase-estimated ZP varied between 0.31 mg C m 23 d 21 and 15.29 mg C m 23 d 21 . The basic multivariate model using biomass, size, and temperature as independent variables, analogous to that presented by Plante and Downing (1989; hereafter PD) for individual populations, performed poorly with our community-level data, explaining only 16% of ZP variation ("model A"; Table 2). While our temperature coefficient (0.05) was roughly similar to theirs (0.03), our body size coefficient was not only different in magnitude (0.05) but also in sign (positive) in our model relative to PD's population-based annual estimate (20.16). Although the general cross-lake trends derived from the PD model agreed with those of our ZP estimates (Fig. 2, r 5 0.53, p 0.001), the PD model systematically overestimated community ZP on average (mean absolute error 5 1.6 mg C m 23 d 21 , an overestimation of 70% on average).
Best subsets regression yielded six models (Table 2, models B-G), and the best of the 1 to 3-variable models were selected using AIC. The best single-variable models (B and C) included temperature and lake area, and were almost equivalent, both describing 12% of ZP variation. The best two-variable models (D and E) included either lake area and temperature, or lake area and Chl a, both explaining 26% of total variation in ZP. The best three-variable models (F and G) incorporated temperature or TP to the previous models, explaining 31-32% of ZP variation. Inclusion of zooplankton standing biomass and mean size with the environmental variables (H) further improved the model fit (AIC 5 17) explaining 37% of ZP variation. Finally, the best predictive model (I) further included % wetland in the catchment (R 2 adj 5 43%, Table 2). The final SEM (Fig. 1b) was significant (p-value 5 0.24, X 2 5 25; df 5 21, CFI 5 0.99) with significant coefficients (all p-values < 0.05, see Supporting Information Table S1) for each relationship. As expected, ZP was positively related to total biomass, and also negatively related to mean body size. The distal (watershed and lake morphometry) and proximal (lake chemistry, temperature) drivers exerted their influence on ZP indirectly, through their effect on zooplankton biomass and body size, and in some cases, also directly, the latter suggesting direct effects on zooplankton performance and growth for any given biomass or body size. The influence of in-lake variables (DOC, TP, and temperature) on ZP was mainly mediated through their effects on Chl a, although temperature also influenced ZP both directly and through an effect on mean body size. The distal drivers, led by % forest in the catchment (%FOR) exerted their influence on ZP through their effect on biomass and mean body size, as well as on in-lake variables, although % wetland cover in catchment (%WET) and lake size were also directly related to ZP and their direct effect was much more important compared to their indirect effects (Table 3). The overall (direct and indirect) effect size of all environmental factors on ZP together was larger (positive 0.98 and negative 20.70; Table  3) than the effect size of community biomass and size alone (positive 0.34 and negative 20.24; Table 3). The SEM also Fig. 2. Comparison of zooplankton production using the PD (Plante and Downing, 1989) model and zooplankton production estimated using chitobiase. The 1 : 1 line is shown on the diagonal. Table 3. Cumulative direct and indirect standardized effects from the SEM explaining ZP using community, within-lake, and catchment-level variables. Variables are: Area, lake area; Biom, standing biomass; Chl a, chlorophyll a; DOC, dissolved organic carbon; %FOR, percent forest cover in catchment; Size, mean individual body mass; Temp, epilimnion temperature; TP, total phosphorus; %WET, percent wetland cover in catchment.

Direct
Indirect 1  demonstrated additional interactions between standing biomass and body size and DOC and TP.

Discussion
There was a surprisingly large range in total daily crustacean ZP estimates across sampled boreal lakes, especially considering that most of lakes were oligo-to mesotrophic. Interestingly, while total biomass also varied greatly across lakes, it was by itself a very poor predictor of daily production, suggesting much variability in biomass-specific production across communities. Previous studies of zooplankton in lakes had reported large differences in daily and annual total crustacean production estimates across individual populations, even within the same lake, with both temperature and body size playing major determinant roles (Banse and Mosher 1980;Downing 1984). Plante and Downing (1989;PD model) further demonstrated that these factors are consistently important across a wide range of aquatic invertebrate species based on empirical relationships, yet no study had tested whether these scalings hold at the communitylevel itself. Our enzymatic approach allowed us to obtain an independent estimate of total daily crustacean production rate that does not rely on population level estimates, and which reflects shorter frames than the latter, more comparable to the scales of variability in the main environmental drivers. Although the enzymatic assay has been itself calibrated using individual, species-specific growth increments, the relationship between enzyme activity and biomass production is consistent across crustacean groups (Sastri and Dower 2009;Sastri et al. 2013), such that the total enzymatic activity (and its turnover) measured in a given water sample offers an integrative quantification of the total daily ZP.
Our results confirm that while community estimates of daily ZP are also dependent on body size and temperature, these dependencies are strongly modulated, and sometimes masked, by several environmental factors appearing to act only at the community level. The direct temperature dependency of ZP was lower than in the PD model, with the direct body size dependency (Model A, Table 2) being positive, and opposite to allometric expectations (Brown et al. 2004). However, after accounting for a suite of environmental factors, including Chl a, lake size, and catchment properties (Model I, Table 2), the resulting body size dependency (M 20.25 ) became negative, as previously found (Stockwell and Johannsson 1997), and of the expected magnitude (Dickie et al. 1987); the temperature dependency was also more in line with previous reports (Plante and Downing 1989;Shuter and Ing 1997) once accounting for additional environmental variables. It would thus appear that both temperature and mean size may integrate the effects of other environmental variables in boreal lakes, such that in order to isolate the temperature and size dependencies of ZP at the community level, it is necessary to account for the effect of these environmental variables. From a mechanistic perspective, these results indicate that body size and temperature dependency of freshwater ZP is modulated by factors (food and other environmental conditions) that likely influence community composition, which then influences the body-size dependency and the response of zooplankton to temperature (Shuter and Ing 1997). The influence of food availability on the body size dependency of ZP had been demonstrated for marine copepods (Lin et al. 2013), but was not considered in previous freshwater, population-based models, which may also explain why the apparent size dependency in the empirical PD model was lower than expected.
Our study identified novel environmental variables affecting ZP at the community level. In their compilation, Plante and Downing (1989) observed no influence of lake area on population-level annual production, while we observed a distinct influence of lake size, with greater ZP in smaller lakes. This discrepancy suggests that lake area acts only at the community level, perhaps influencing community composition or species interactions along a gradient of lake morphometry. Smaller lakes are more likely to be fishless (Scheffer and Van Nes 2007), resulting in greater invertebrate vs. fish planktivory, with important impacts on zooplankton community composition (Brooks and Dodson 1965;Hanazato and Yasuno 1989). Compositional alteration can in turn lead to life history shifts affecting reproductive rates that modulate community ZP without necessarily being accompanied by altered mean community body size. We also found that ZP was greater in catchments with lower % wetland coverage. While we were unable to identify the intermediary inlake variable underlying this relationship, it is well known that the presence of wetlands in the catchment strongly influences material loading to lakes, notably N, P and especially DOC (Dillon et al. 1991), influencing water chemistry and concomitantly, the composition of phytoplankton (Prepas et al. 2001). It is interesting to note that the wetland effect in the multiple regression models was not diminished by the inclusion of variables such as nutrients and total or colored DOC, indicating that this % wetland likely integrates several different processes influencing zooplankton community structure, not just water color.
The SEM approach revealed more of the interactions of environmental and biological variables shaping community ZP. As expected, the effect of within-lake variables on ZP was mainly mediated through food concentration (Chl a), although direct and indirect (through body size) temperature effects were apparent. Catchment % forest cover, a watershed level variable similar to % wetland, indirectly influenced ZP via both in-lake variables and community properties, likely represents the integration of several watershed processes influencing lake chemistry, resource availability and even physical structure that affect zooplankton community composition and interactions. Although the causal link to ZP is not clear, it is interesting that the standardized cumulative effect of catchment-level variables F. St-Gelais et al.
Zooplankton community production was similar to the cumulative effect of within-lake variables, emphasizing a role of catchment properties on regulating daily zooplankton production. Similarly, while lake area and temperature mainly exerted direct effects on ZP, indirect effects were also apparent, operating through mean body size. Together, these results imply a complex regulation of daily zooplankton production at the community-level. While the Plante and Downing (1989) population-based model applied to our community data yielded estimates that were greater than our measured ZP, the resulting patterns were still generally consistent with our community-level ZP observations (Fig. 2). This would suggest that that the basic scaling of zooplankton production to biomass, body size, and temperature is roughly comparable between the population and community levels. Our results further demonstrate, however, that community estimates capture an additional layer of cross-lake variability, which is at least in part explained by a set of integrative watershed and lake drivers that do not appear to operate at the population level. Our results clearly indicate that total daily crustacean production regulation should not be assumed to be simply the sum of its constituents, and that regulation occurs through a complex interplay between community structure, lake, and watershed properties.