Environmental drivers of mixotrophs in boreal lakes

Mixotrophy is increasingly recognized as an important trophic pathway among phytoplankton, yet its underlying drivers remain largely unknown and unexplored. Here, we present a study utilizing 69 lake samplings in boreal Quebec, Canada, identifying variables driving the success of phytoplankton that have a capacity for mixotrophy and pointing to the underlying mechanisms. We found that the success of mixotrophs (% of total biomass) was positively influenced by both colored dissolved organic matter (cDOM) and dissolved CO2 concentration but limited by the abundance of crustacean zooplankton. The effect of cDOM manifested as a consequence of limited autotrophic phytoplankton biomass in lakes with reduced light penetration. We observed a nonlinear (u‐shaped) relation between CO2 and mixotrophs, with biomass favored at both low and high CO2 concentrations. A reduced fitness of mixotrophs at near‐atmospheric CO2 concentrations is likely owing to the costs of rapidly switching between or maintaining multiple trophic strategies. The abundance of zooplankton had a negative effect on mixotroph biomass but a positive effect on autotrophic phytoplankton. We also found that while the community composition of potentially mixotrophic phytoplankton was to some degree likely influenced by zooplankton biomass, composition was unaffected by the CO2 and cDOM gradients. Overall, this study highlights mixotrophy in boreal lake systems as a strategy of persistence, with the maintenance of a moderate but constant presence across a changing gradient of light and trophic conditions. The results of our study support the hypothesis that phytoplankton with a capacity for mixotrophy provide a superior and stable stoichiometric food source for zooplankton, implicating mixotrophs as a vital component of boreal lake food webs.

Mixotrophy among aquatic microorganisms such as planktonic protists is most frequently defined as the ability to obtain energy and nutrients through a combination of photosynthesis and heterotrophy. This classifies mixotrophy as a generalist type of trophic strategy (Selosse et al. 2017). As a functional trait, a capacity for mixotrophy has long been viewed as a fringe strategy in lake and ocean ecosystems. However, recently, mounting evidence has shown this to be a severe underestimation with examples of phagotrophy found in most phytoplankton groups, indicating that mixotrophy plays a major role in lake carbon cycling and food web dynamics (Flynn et al. 2013;Mitra et al. 2014), with studies indicating that some mixotrophs may account for more than 60% of bacterivory in certain aquatic ecosystems (Unrein et al. 2007;Zubkov and Tarran 2008).
Phytoplankton with the capacity to obtain energy and critical nutrients through a combination of innate photosynthesis (not acquired through ingestion of phototrophic prey) and phagotrophic bacterial grazing have been defined by Mitra et al. (2016) as constitutive mixotrophic phytoplankton.
An important distinction to make in the study of mixotrophy is between estimation of rates of nutrient and energy acquisition through multiple trophic pathways (i.e., the act of mixotrophy) and that of the capacity for mixotrophy as a trait (i.e., the ability to switch between trophic pathways or utilize multiple pathways simultaneously, to varying degrees). Some studies utilize the former definition, considering mixotrophy rates in groups of organisms at the moment and location where a study occurs (Ptacnik et al. 2004(Ptacnik et al. , 2016Fischer et al. 2017). However, this view excludes groups of organisms that may not be actively utilizing their mixotrophic capacity at the precise moment of study. It also overlooks the importance of mixotrophy as a plastic trait that may enable an organism or a group of organisms to persist under a wider range of environmental conditions. The approach we took therefore relied on relating the group-specific (pure autotroph or mixotroph) biomass of phytoplankton to environmental factors (Reynolds et al. 2002).
Combining different trophic pathways has some intuitive advantages, with the most obvious one being that mixotrophs may substitute limited resources with alternative ones, or to obtain them through a more efficient pathway according to the environment (Bird et al. 1998;Isaksson et al. 1999). The resources acquired vary between types and species of mixotrophs, and mixotrophy can be a way of fulfilling either nutrient (phosphorus/nitrogen) or carbon requirements (Stoecker 1998), whereas for others it is a means to obtain specific metabolites, for example, phospholipids (Kimura and Ishida 1989). However, the costs of maintaining an ability to switch between trophic strategies or to perform them simultaneously are bound to be significantly higher than maintaining only one specialized feeding pathway at all times (Raven 1997;Flynn and Mitra 2009;Andersen et al. 2015). Thus, mixotrophs should be at a disadvantage relative to heterotrophs either under low light conditions or at high bacterial prey availability (Fischer et al. 2017). Conversely, they should also be disadvantaged relative to photoautotrophs under low light conditions, when nutrients are replete, as photoautotrophs compete better under low light owing to their higher cellular chlorophyll contents (Rothhaupt 1997;Tittel et al. 2003;Ptacnik et al. 2004Ptacnik et al. , 2016Katechakis and Stibor 2006;Bruggeman 2011;Fischer et al. 2017). In contrast, mixotrophs should generally be favored in environments where either specific resource requirements (e.g., nutrients) cannot be met by photoautotrophs or prey concentration for obligate phagotrophs is too low to cover their energy and nutrient demands while foraging for prey (Ptacnik et al. 2016;Fischer et al. 2017). In the latter scenario, mixotrophs have an advantage by being able to sustain themselves through added energy from photosynthesis (Fischer et al. 2017). Such conditions exist in clear, oligotrophic lakes, where light availability provides energy to mixotrophs that supplements the extra energy costs of maintenance, and bacterivory supplements the limited availability of nutrients that is detrimental to pure photoautotrophs (Calbet et al. 2012;Ptacnik et al. 2016). As such, mixotrophic success is currently, to a large part, understood to be related to a specific set of environmental conditions: low nutrient, but high light availability.
Such environmental limitations stand in contrast with recent findings indicating that mixotrophs are much more commonly distributed than previously thought, occurring over a range of environmental conditions (Flynn et al. 2013). Their widespread distribution is also in contrast to the idea that maintaining a mixotrophic potential carries disadvantages in the form of added structural and physiological costs. Although mixotrophs should not be able to directly compete with either photoautotrophs or heterotrophs under conditions of high nutrient or low light availability, they may be superior when multiple or all resources are scarce by pooling the energy and nutrient uptake from multiple sources, thereby overcoming limitations by the relatively low abundances of each resource and the associated costs. Indeed, several other potential drivers of mixotrophy have been revealed recently. For example, Bergström et al. (1998) found that mixotrophs dominated the phytoplankton communities of three Swedish brown lakes, suggesting that light limitation through humic matter from the catchment area could potentially favor mixotrophy, an observation confirmed by several other studies of similar lakes (Ramberg 1979;Rothhaupt 1996a,b;Bergström et al. 1998;Bergström et al. 2000;Bergström 2009). In a study where thermocline depth was experimentally altered in a lake, thus leading to hypolimnetic (limiting) nutrient access in darker waters, mixotrophic taxa, especially cryptophytes, were the noted beneficiaries (Cantin et al. 2011). These empirical studies indicate that while the prevalence and relative success of mixotrophs in light-rich, nutrient-limited environments (i.e., clear water and oligotrophic lakes) have been well-documented, this may only provide a glimpse into a larger picture; one in which mixotrophs may act as opportunists or generalists, rather than specialists. Taken together, these studies suggest a much wider ecological role of mixotrophs as a group, that prevails in a wide range of habitats by maintaining a more stable and constant presence than their more specialized competitors, allowing them to benefit from opportunities arising from punctuated environmental shifts such as humic input or water column mixing.
Interestingly, while there are several studies addressing the abiotic drivers of mixotrophy such as light and nutrients, few have investigated other potential drivers, such as predation, and how such drivers may influence the outcome of their competition with other phytoplankton. Indeed, owing to their flexibility in nutrient acquisition (Gächter and Bloesch 1985;Elser and Hassett 1994;Elser et al. 2000;Katechakis et al. 2005), it has been demonstrated that mixotrophs at least under nutrient-poor conditions may have a higher or more stable nutritional value than autotrophic phytoplankton. As zooplankton tend to be highly affected by stochiometric composition of their prey (Hessen 1992;Gulati and Demott 1997;Hessen et al. 2003), it has been hypothesized that zooplankton may prefer mixotrophic phytoplankton, at least under certain environmental conditions (Ward and Follows 2016;Katechakis et al., 2005). To this end, Katechakis et al. (2005) in a chemostat experiment found a preference among zooplankton for at least certain mixotrophs, whereas others appeared to be toxic.
To date, although studies have indicated dominance of mixotrophs under light-limited conditions, there has been limited research addressing the prevalence of mixotrophy compared to photoautotrophy under conditions other than those already confirmed to favor a mixotrophic strategy (i.e., limited nutrient and high light availability; e.g., Calbet et al. 2012;Ptacnik et al. 2016) or indeed, comparing the relative success of the phytoplankton trophic groups in general (Jones 2000). In this study, we aim to investigate whether mixotrophic success in situ can be driven by other environmental variables, including light limitation via browning, i.e., high concentrations of colored dissolved organic matter (cDOM), or zooplankton predation. We expect that mixotrophic success in situ may also be driven by these factors in several ways. First, we hypothesize that greater light attenuation caused by cDOM will limit photoautotrophic phytoplankton growth to a larger extent than mixotrophic growth, thereby leading to an increased proportional success of mixotrophs (Hypothesis 1) as they may maintain biomass by switching to more heterotrophic feeding. Second, because zooplankton are thought to benefit from consuming mixotrophic prey (based on nutrient stoichiometry; Katechakis et al. 2005), we hypothesize that evidence for top-down control by zooplankton on mixotrophs and their contribution to overall phytoplankton biomass should be present (Hypothesis 2). To test these hypotheses, we examined phytoplankton communities, considering the absolute and relative biomasses of potentially mixotrophic and known photoautotrophic taxa, in relation to several limnological factors across 55 North American lakes.

Study sites
A total of 69 samples collected from a set of 55 lakes with a low level of anthropogenic disturbance were sampled within three environmentally and geologically distinct regions of Quebec, Canada (Fig. 1). Although some lakes were sampled twice, this was done at different locations, at an interval of 2 months or more. Because of the difficult access to such remote lakes, sampling occurred in three different summers for each region: Abitibi in 2010, Chicoutimi in 2011, and Schefferville in 2012. For overview of variables and variable explanation, see Table 1.

Biological variables
The phytoplankton protist communities were sampled at the deepest point in each lake over the photic zone using a flexible poly(vinyl chloride) sampler tube. An integrated subsample (250 mL) was preserved in Lugol's solution, and phytoplankton were enumerated in the lab at the genus level using the Utermöhl method on an inverted microscope at 400X magnification until a minimum of 300 individuals or colonies were counted in separate fields of view. Unfortunately, heterotrophic nanoflagellates and ciliates could not be identified from these samples. Phytoplankton biomass by taxon was estimated from biovolume computed using cell and colony length measurements and corresponding geometric shapes (Hillebrand et al. 1999). Mixotrophic phytoplankton taxa were assigned in our dataset based on a literature review (Longhi and Beisner 2010), focusing on known photosynthetic bacterivores (primarily cryptophytes, chrysophytes, dinophytes; chlorophytes were considered autotrophic). Estimates of the total biomass of mixotrophs (Mixo_Biom) as well as their proportion of total phytoplankton biomass (%Mixo) were calculated. The total biomass of purely autotrophic phytoplankton was also estimated (Auto_Biom).
To obtain total zooplankton biomass, integrated samples were also collected from the deepest point in each lake using a conical plankton net (110 μm mesh, 0.30 m mouth diameter), equipped with a flow meter (General Oceanics), hauled vertically from 1 m above the sediments to the surface. Zooplankton samples were anesthetized using carbonated water and were preserved in 75% (final concentration) ethanol. Crustacean zooplankton species were identified using an inverted microscope (50-400X) and enumerated until a total of 200 individuals had been identified. For each taxon present in each lake, the length of 20 mature individuals was measured and biomass by taxon was estimated using length-dry-mass regressions (McCauley 1984;Culver et al. 1985).

Limnological data
Key limnological variables were measured to characterize the lake and catchment environments. We used a multiparameter sonde (Yellow Springs Instruments [YSI]) to measure pH (at 0.5 m) and temperature (at 0.5 m, then averaged over the water column). Water samples were collected at 0.5 m to measure the concentration of chlorophyll a (Chl a), total phosphorus (TP), total nitrogen (TN), bacterial abundance, dissolved organic carbon (DOC), cDOM (an indicator of lake color), and fluorescent DOM (fDOM). Chl a was extracted with 90% hot ethanol and absorption was measured spectrophotometrically before and after acidification to account for phaeophytin after (Lorenzen 1967;Mush 1980); TP was measured from water samples using the molybdenum-blue method following persulfate digestion (Cattaneo and Prairie 1995); TN was measured using nitrates after persulfate digestion; bacterial abundance was measured by cytometry (FACSCalibur, Becton Dickinson). DOC concentrations were measured on an OI 1010 TIC-TOC analyzer following sodium persulfate digestion, and cDOM was quantified as the absorbance at 440 nm using an Ultrospec 3100 spectrophotometer (Biochrom) followed by PARAllel FACtor analysis (PARAFAC) to decompose the fluorescence signals into their underlying chemical components (Murphy et al. 2010(Murphy et al. , 2013. In the ideal case where fluorescence conforms to Lambert-Beers law, this process can lead to the mathematical identification and quantification of independently varying fluorophores (Murphy et al. 2010(Murphy et al. , 2013Stubbins et al. 2014). We used PARAFAC-identified components as a measure for different DOM types. Components 1-4 (fDOM_1 to fDOM_4) represent humic/fluvic acids, typically of terrestrial origin, whereas components 5-6 (fDOM_5 to fDOM_6) are associated with more protein-like molecules associated to a high biological degradability (Stubbins et al. 2014). Lake area was derived using ArcGIS V10 software (ESRI). CH 4 was measured as the mean methane concentration in water (parts per million), extracted by headspace (3 × 60 mL syringes) and measured by gas chromatography (GC-2014). Dissolved CO 2 concentration was measured as mean carbon dioxide concentration in water (parts per million), extracted by headspace (3 × 60 mL syringes) and measured with EGM-4 (PP Systems). Lake respiration was calculated as mean community respiration (μg C L −1 d −1 ) measured by O 2 depletion in incubations (Fibox, PreSens). Latitude and longitude of the sample site (decimal degrees) measured by Global Positioning System (GPS) and altitude were also determined, maximum depth measured at each sampling site with a sounding line. Mixing depth was measured as depth of epilimnion (meters) determined from vertical profiles of temperature and O 2 , measured with YSI, and CO 2 measured with an EGM-4. Secchi depth (m) was also measured at each sampling site.
When there was a missing value for the limnological data (occurred once for pH, once for water temperature, DO, and TP, twice for Chl a and Zoo_Biom, thrice for maximum depth, and four times for thermocline depth), values were imputed using a Random Forest (RF) approach (missForest R package Stekhoven and Bühlmann 2012, NRMSE: 0.038). To avoid including highly correlated explanatory variables, prior to further analyses, we first excluded those with r >AE 0.7, removing the variable that was least correlated with the %Mixo response variable. Variables excluded in this way were cDOM, due to high correlation with fDOM_4, as well as DOC, fDOM_2, and fDOM_3, due to high correlation with fDOM_1. Because fDOM_4 and cDOM (r = −0.66 vs. r = −0.65, p < 0.001 in both) correlated similarly with Secchi depth, we used fDOM_4 as the most direct indicator of DOM-induced light limitation. Using PARAFAC measurements as a proxy for DOMinduced light limitation for phytoplankton by specific DOM compounds is to our knowledge a novel approach, which may yield more specific information with regards to properties and effects of DOM on these communities.

RF analyses
We used RF analysis to analyze the influence of preselected variables in the larger dataset on the three main response variables (Mixo_Biom, Auto_Biom, and %Mixo). This was also used in a secondary set of analyses to reveal potential drivers of CO 2 , a nonlinear explanatory variable of phytoplankton. RF has several advantages over traditional regression techniques. It is considerably less vulnerable to overfitting than is multiple regression when processing a large number of predictor variables. Furthermore, RF does not assume a specific distribution in the data, is highly robust toward outliers, and can identify effects of advanced interactions, nonlinear relationships, and even nonmonotone association rules, without these being specified beforehand (Matsuki et al. 2016;Ryo and Rillig 2017).
The RF algorithm is a machine-learning technique based on an ensemble of regression trees (Breiman 2001). RF bases each tree on a randomly selected subset of the total number of variables, before it chooses the most influential variables based on consensus from its regression tree ensemble. A response class is predicted in each terminal node of the tree (or each rectangular section in the partition, respectively) by deriving from all observations in this node either the average response value in regression or the most frequent response class across all classification trees. Note that this means that a regression tree creates a piecewise (or rectangle-wise for two dimensions and cuboid-wise in higher dimensions) constant prediction function. This decreases the correlation between trees, and therefore also reduces the error rate of the RF as a whole. Here, we used a RF technique based on conditional inference regression trees (Strobl et al. 2009a) developed by Ryo and Rillig (2017). A measure of importance was calculated for each predictor variable by cross validating each tree with data not used when the tree was constructed, referred to as the out-of-bag (OOB) data (Breiman 2001).
We conducted separate RF analyses for each of the three main response variables. In each analysis, 5000 regression trees were used to obtain a stable prediction. A large number of trees was necessary, because despite removing highly correlated variables (r >0.7), there were still several explanatory variables that had high correlations (r >0.6), for example, pH and CO 2 . The RF from party package was chosen specifically, because it is better at dealing with such issues (Hothorn et al. 2006;Zeileis et al. 2008;Strobl et al. 2007;Ryo and Rillig 2017). However, while a general division between more and less important variables could be established with a smaller number of forests (500), a moderate number of >5000 regression trees was required to obtain a relatively stable tree.
To simplify categorization of variable importance, and because certain variables were considerably more influential than others, we defined those with major influence as explaining more than 10% of the portion of variation explained by the RF. Variables explaining 5-10% were defined as being of intermediate and variables of 1-5% as of low influence (<1 are considered of minor influence, and were considered too small to be reliable due to inherent instabilities in the RF from run to run). Due to potential accumulation of uncertainty in each analysis, we refrained from discussing variables which did not exceed 10% in at least one of the three RFs.
To visualize our results, we created partial dependency plots (PDPs) of the relationships between each of response variable against scoring predictor variable from the RF. PDPs do not display the data directly but are projections based on the model inferred by the RF. As such, the range of the variation displayed on the y-axis is the proportion of the range explained by the variable in question (according to the RF) (Carlisle et al. 2009;Leach et al. 2018). Relations between community and environmental variables were also explored using regressions. Independent regression plots from multiple response variables (Mixo_Biom + Auto_Biom + %Mixo) were superimposed against each driving variable to examine directional variation between these phytoplankton community responses. Because linear regression (unlike RF) is sensitive to outliers, nonlinearity and non-normal distributions, variables with such issues were transformed, rather than having outliers removed.

Ordinations of community composition
To determine the degree to which changes in community composition drives the output of our RF analyses, we conducted a series of nonmetric multidimensional scaling (NMDS) and redundancy analysis (RDA) using the biomasses of dominant phytoplankton genera (genera occurring at more than three sites) and zooplankton as response variables. For genera of mixotrophs and autotrophs, see Table 2. Zooplankton were classified as: Daphnia, other Cladocera (OtherClads), Cyclopoid Copepods (CycCopes), and Calanoid Copepods (CalCopes). NMDS analyses using the metaMDS function from the vegan package in R (2018; R package vegan) were conducted separately for biomass of mixotrophs, autotrophs, and zooplankton across lakes, according to a classification of lakes as being either of low or high concentration of CO 2 (low <425 ppm and high >1040). Intermediate CO 2 concentration lakes were excluded. Similarly, NMDS was done by classifying lakes as of either low or high Zoop_Biom type (low <41,570.56 μg m −3 and high >197,136.92 μg m −3 ) and low or high fDOM_4 type (low <0.1 and high >0.193). The high and low concentration groups were defined based on the PDPs for Mixo_Biom against each of these explanatory variables. To test for significance, the NMDS was followed by a permutational multivariate analysis of variance (PERMANOVA) test using the adonis function in the R package vegan.
Separate RDAs (using the RDA function from vegan package) were done to evaluate variation in the biomass composition of mixotrophic and autotrophic phytoplankton and zooplankton, with the set of explanatory variables based on results from the RF for %Mixo (Table 3). To compensate for non-normal distributions, and strong outliers (defined as data points beyond two quantiles from the median), CO 2 , pH, and Zoop_Biom were log transformed. Biomass matrices were Hellinger transformed before analyses. Permutation (999) tests were used to test the significance of each model, axis, and terms in each RDA using the anova function of the vegan package in R. To examine potential relationships between mixotrophs and zooplankton predators, we also constructed an RDA wherein Zoop_Biom in the explanatory variable data matrix was replaced with the main zooplankton group biomasses (Daphnia, OtherClads, CycCopes, and CalCopes).
All data analyses were carried out in the statistical programming environment R, version 3.4.2 (R Development Core Team 2017). The RF analyses were done with an implementation written by Strobl et al. (2009b) (R package party).

Interpretation
Because RF is a more robust and reliable method than linear modeling when it comes to measuring influential relationships, we were able to use the results of the RF to inform other linear analyses (simple regression, RDA and NMDS) to draw conclusions. RF provided ranked information on influential variables of phytoplankton community composition. We then used regression plots in combination with the PDPs (RF) to discern the directional relationships and compare these between response variables. RDA and NMDS were used to examine responses and effects at the genus level and to determine whether observed effects also represented shifts in community composition within the mixotroph and autotroph groups.

Random forests
The RF analyses including all the selected variables explained 51% of the variation in %Mixo for all the sampled lakes, with an OOB error rate of $1% in the test dataset. The model explained 42% of variation in the mixotroph biomass (Mixo_Biom; OOB 2%) and 37% of variation in the autotroph biomass (Auto_Biom; OOB 3%) ( Table 3).

Relationship directionality
Examination of the PDP (Fig. 2) and linear regression plots (Fig. 3) provided insight into the direction of the relationships of %Mixo with the factors identified in Table 3. %Mixo was strongly and positively related to both fDOM_4 and Slope, occurring via a strong negative effect on Auto_Biom with no detectable effect on Mixo_Biom. Conversely, %Mixo declined  Table 1.
with Zoo_Biom via small negative effects on Mixo_Biom but positive on Auto_Biom. %Mixo was strongly and positively related to dissolved CO 2 , occurring via strong positive responses of Mixo_Biom and a strong negative one of Auto_Biom. There was a negative relationship between %Mixo and pH, driven by a stronger positive effect on Auto_Biom than on Mixo_Biom. Finally, %Mixo declined with Bact_Resp, via smaller relative decreases in the Auto_Biom portion of the phytoplankton communities than in the Mixo_Biom portion. No net effects on % Mixo were observed for Chl a, Water_Temp, and Bact_Abund, which all positively influenced both Auto_Biom and Mixo_ Biom. However, Bact_Abund was somewhat more related to Mixo_Biom.
Factors explaining CO 2 RF on CO 2 as response variable gave three potential main drivers in order of importance: fDOM_4 (variable importance, VI = 24.0), TP (VI = 23.3), and Secchi (VI = 10.7) (Supporting Information Fig. S1). Below CO 2 = 600 ppm, no relationships were determined in the RF, which might explain the relatively high OOB = 21.1% despite an overall model R 2 of 57%. At higher CO 2 concentrations, CO 2 increased nonlinearly with TP and fDOM_4 and decreased with secchi and DO.

NMDS of community composition under extreme driver values
NMDS followed by PERMANOVA revealed no significant differences between mixotroph biomass community composition at extreme (low vs. high) CO 2 values (Supporting Information Fig. S2). Conversely, autotrophic phytoplankton compositions did differ (p = 0.018), although only a small fraction of variation (R 2 = 0.08) was explained by the model (Fig. 7A) (C) 6.5 7.0 7.5 8.0
There were no significant differences in zooplankton community composition at extreme values (low vs. high) of CO 2 (Supporting Information Fig. S6) or fDOM_4 (Supporting Information Fig. S7). However, zooplankton composition did vary between the two Zoo_Biom levels (p = 0.003, R 2 = 0. 225), with daphniids dominating lakes that had high total zooplankton biomass (Fig. 8).

Discussion
Our study of 55 North American lakes revealed that mixotrophy was greatly controlled by environmental variables such as cDOM ("fDOM_4"; Table 1), which influenced light availability, as well as food web composition including heterotrophic  activity and crustacean zooplankton. Our results will help to inform future studies about the occurrence and ecological role that mixotrophy plays in aquatic food webs under varying regressed against %Mixo RF variables. x-axis = RDA1, y-axis = RDA2. Constrained and total variation explained by RDA1-axis and RDA2-axis is given in each axis legend, respectively. For an overview of variable units and variable description, see Table 1. Zooplankton were classified as: Daphnia, other Cladocera (OtherClads), Cyclopoid Copepods (CycCopes), and Calanoid Copepods (CalCopes). [Color figure can be viewed at wileyonlinelibrary.com] environmental conditions. The overall success of mixotrophs, defined as proportion of the phytoplankton community composed of potentially mixotrophic taxa, was mainly linked to increases in fDOM_4 and dissolved CO 2 concentrations and reduced zooplankton biomass. Positive effects of Slope (i.e., average downward angle of the catchment area) and negative effects of pH were of more intermediate importance for this indicator of mixotrophic success. Conversely, while water temperature and bacterial abundance both positively influenced the standing biomass of mixotrophs in lakes, they did so in the same relative positive way for the autotrophic phytoplankton. Overall, our results indicate that the proportion of mixotrophy is higher in lakes receiving more allochthonous DOM (highly correlated with the fDOM_4 indicator; r = 0.84, p < 0.001), but where zooplankton predation pressure is low. Taken together, the importance of mixotrophy in brown water and oligotrophic lakes manifests as a result of interactions between abiotic and biotic drivers, implying a better understanding, can be gained from observational studies.

Lake cDOM and browning favor mixotrophy
A recurring problem when dealing with cDOM is that it is a parameter which may consist of a myriad of different compounds from different sources and with different properties (Stubbins et al. 2014(Stubbins et al. , 2017Wagner et al. 2015). By using a PARAFAC measurement rather than cDOM directly, we are able to more accurately identify the DOM components influencing mixotrophic and autotrophic phytoplankton (in the case of fDOM_4, fresh, allochthonous, and labile substances). This study is to the best of our knowledge, the first using PARAFAC as a proxy for lake cDOM to study effects on phytoplankton. Influences of fDOM_4, Slope, and CO 2 likely represent a suite of factors that favor mixotrophy, most likely via numerous linked processes related to allochthonous DOM input into lakes (Bergström 2009). Because fDOM_4 was the PARAFAC component that most strongly correlated with cDOM and matched its correlation with Secchi depth, we interpret this as fDOM_4 being the cDOM component responsible for light limitation (as increased cDOM does not always equate to decreased light penetration). Highly labile fDOM_4 (as is the case for the fraction selected by the RF model: fDOM_4; Stubbins et al. 2014) possibly augmented bacterial biomass, thereby increasing CO 2 levels in the water column.
No correlation was observed between fDOM_4 and Bact_Abund; as would be expected given that an increase in the topdown mixotrophic grazing pressure would mask a positive bacterial response to fDOM_4. The positive relationships between Bact_Abund and both Mixo_Biom and Auto_Biom are most likely caused by a positive feedback on the bacterial community via photoautotrophic excretion of labile products, while Bact_Abund exerted a positive effect on mixotrophs as a potent food source (Fig. 3;Bauerfeind 1985;Simon and Tilzer 1987). While all of these variables were associated to greater % Mixo (Bact_Abund excluded), notably they did not directly support greater mixotroph biomass, as at the same time light availability for autotrophs was reduced (Table 3; Fig. 3). The net effect, however, was an increase in the proportion of mixotrophs, indicating greater competitive success. Thus, it appears that while photoautotrophic activity suffers in "browner" lakes, phytoplankton that adopted a mixed feeding strategy including heterotrophy do not suffer much from reduced light availability (increased fDOM_4). Most likely, mixotrophs engage heterotrophy to access alternative energy sources under reduced light conditions. This is in line with our first hypothesis that reduction in light availability, e.g., with lake browning, will reduce growth of photoautotrophs more than mixotrophs, which maintain a constant presence across a wide range of environmental settings. These findings also correspond with previous studies indicating that mixotrophic phytoplankton is more abundant and even dominates in brown, humic matter-rich, and lightlimited lakes (Rothhaupt 1996a,b;Bergström et al. 1998;Bergström et al. 2000;Bergström 2009). It has been argued that bacterivory allows mixotrophs to outcompete photoautotrophs under light-limited conditions (Jones 2000). A notion which is supported by Jansson et al. (1996) who found strong hints that mixotrophs in humic lake enclosures in Sweden were indeed limited by bacteria, which also are dependent on allochthonous DOC and thus additionally comprise a light limiting component. This indicates a complex interdependent relationship (Bergström 2009). Roberts and Laybourn-Parry (1999) found that mixotrophs were able to survive under the ice cover of Antarctic lakes even in the absence of light in winter, again likely by drawing on alternative energy sources (i.e., bacterivory), while pure autotrophs became limited by one or two primary resources, a finding which has later been confirmed in several studies (Laybourn-Parry 2002;Bielewicz et al. 2011;Dolhi et al. 2012). To our surprise, only a few studies have been conducted to specifically address this hypothesis directly. The majority of literature addresses light availability in combination with nutrient limitation, concluding that under high light and low nutrient availability, mixotrophic protists outcompete photoautotrophs that suffer by a lack of dissolved nutrients (Calbet et al. 2012;Ptacnik et al. 2016;Fischer et al. 2017). To our knowledge, the present study is one of the first to examine, albeit via a low temporal resolution survey of multiple boreal lakes, for the simultaneous effect of light limitation and higher bacterivory along the cDOM gradient.
CO 2 as a suitable indicator of the degree of mixotrophy CO 2 had a significantly positive effect on %Mixo and was one of our major influential variables. As with fDOM_4, the relationship appears to arise from a stronger negative effect on Auto_Biom than on Mixo_Biom. However, the relationship with Mixo_Biom (Fig. 2B) points to a more complicated relationship. Closer examination revealed that both high and low levels of CO 2 appear to favor high levels of Mixo_Biom, whereas at intermediate levels (concentrations close to atmospheric equilibrium), Mixo_Biom was lower. Vogt et al. (2017) examining these same lakes found that gross primary production showed the same pattern, being the lowest at atmospheric equilibrium relative to undersaturated or supersaturated CO 2 conditions. Previously, it has been demonstrated that CO 2 in boreal lakes is closely related to heterotrophic activity, increasing with cDOM (Jonsson et al. 2003) that should induce light limitation, thereby reducing autotrophy. In line with this study, we found a strong relationship between CO 2 and fDOM_4 (Supporting Information Fig. S1). This is probably because CO 2 initially is driven mainly by the degree of autotrophy, which in turn is limited by fDOM_4. As autotrophy declines, CO 2 increases, until autotrophs reach a lower threshold. Mixotrophs appear to follow a similar relationship. However, as CO 2 concentrations attain levels approximately at equilibrium with the atmosphere, mixotroph biomass starts increasing again, eventually reaching values closer to those observed at low (<400 ppm) CO 2 concentrations. Variation in CO 2 concentration was best explained by fDOM_4 (Supporting Information Fig. S1). This could mean that even as autotrophy becomes light limited, increasing fDOM_4 and TP may continue to fuel an increasing concentration of CO 2 (Supporting Information Fig. S1). The labile properties of fDOM_4 coupled with higher TP could fuel lake respiration through heterotrophic production, increasing CO 2 production indirectly. Unfortunately, without data on heterotrophic plankton, we are unable to investigate this in any more meaningful detail. However, in a marine mesocosm study by de Kluijver et al. (2013), mixotrophic phytoplankton showed a slight positive relationship with increasing concentrations of CO 2 in combination with nutrient addition. De Kluijver et al. (2013) attributed this effect to reduced competition with other phytoplankton, processes that could also explain increases in mixotroph biomass that we observed following nearatmospheric CO 2 concentrations.
Whatever the underlying mechanisms of lake CO 2 concentration, together with our observations, this variable emerges as a good proxy for lake primary productivity and the degree of mixotrophy in these boreal lakes. This suggests mechanistically that there is a cost associated with being simultaneously autotrophic and heterotrophic or for rapidly switching between these trophic modes as has been suggested by others (Raven 1997;Flynn and Mitra 2009;Andersen et al. 2015). This is manifested by a reduced ability of mixotrophs to survive and/or compete with purely autotrophs when light availability is high (undersaturated CO 2 lakes). As conditions become more heterotrophic (supersaturated CO 2 lakes), mixotrophs are able to switch to conducting heterotrophy and outcompete autotrophic specialists under light-limited conditions.
Overall, it appears that the success of the mixotrophic strategy in brown water lakes (high fDOM_4) is linked to their resilience to changing conditions as manifested by the maintenance of moderate abundance under a wide range of circumstances, rather than a maximum biomass when environmental conditions allow for it. This is also supported by the fact that we observed no significant variation in mixotroph community composition between low and high CO 2 concentration lakes, indicating that the same mixotroph taxa are able to persist over a range of environmental conditions. This notion further supports our first hypothesis, providing indirect evidence that mixotrophs may sustain their biomass levels (but not increase) against pure heterotrophs (unmeasured in our study) even in the more brown, cDOM-rich lakes.

Zooplankton predation may control mixotrophic success
The relative success of mixotrophy declined with zooplankton biomass, suggesting a strong top-down effect. The net negative effect was a result of opposing responses by the mixotrophic and autotrophic portions of the phytoplankton community to zooplankton: Auto_Biom responding surprisingly positively, while Mixo_Biom responded negatively. These results point to a preference of the zooplankton community for mixotrophic over autotrophic phytoplankton prey, supporting our second hypothesis that zooplankton should limit the contribution of species with the propensity for mixotrophy to overall phytoplankton biomass.
Evidence for zooplankton preference of mixotrophic phytoplankton is limited to date, and studies addressing this subject are somewhat conflicting, with preference for certain species of mixotrophs confirmed while evidence of toxic effects are shown in others (Katechakis et al. 2005). Advanced prey selection by zooplankton appears to occur via chemosensing, as well as other senses based on phytoplankton morphology and behavior (DeMott, W.R., 1986;Broglio et al. 2004;Katechakis et al. 2005). However, a general preference for mixotrophs by zooplankton could be related to elemental nutrient stoichiometric ratios, while toxic effects may be a result of coevolution resulting in a predator defense. Elemental composition is a key factor influencing food web trophic efficiency, because there are often mismatches between carbon : nutrient ratios of consumers vs. producers (Hessen 1992;Gulati and Demott 1997;Sterner and Schulz 1998;Hessen et al. 2003). In particular, C : element ratios in autotrophs can often be suboptimal for zooplankton grazers as C : P ratios of many photoautotrophs may be too high and too variable with respect to the requirements of herbivorous zooplankton (Gächter and Bloesch 1985;Elser and Hassett 1994;Elser et al. 2000). Furthermore, the light nutrient hypothesis purports that biomass and nutrient stoichiometry of zooplankton that feed primarily on photoautotrophs in particular depend on light and P-supply (Andersen 1997;Sterner et al. 1997;Loladze et al. 2000) such that herbivore growth and fecundity are limited by food quantity at low light intensities and by stoichiometric food quality at high light intensities. Conversely, nutrient stoichiometry of mixotrophs has been found to remain almost constant even when ratios of light : dissolved nutrients vary (Gächter and Bloesch 1985;Elser and Hassett 1994;Elser et al. 2000;Katechakis et al. 2005). Additional studies have found that zooplankton benefit from a stable intake of polyunsaturated fatty acids in their diet (Brett and Müller-Navarra 1997;Ravet et al. 2003), while other separate studies have found some indications that certain mixotrophs may be able to maintain relatively stable and/or high levels of such substances (Garcí et al. 2000;Adolf et al. 2007). Thus, mixotrophs may supply zooplankton with a more stable diet, especially when light : nutrient regimes are variable, which provides a mechanism for the apparent preference of zooplankton for mixotrophs in our study.
Analyses of community compositions through RDA revealed that mixotrophic species composition, to some degree, is driven by the slope of the landscape (with higher values of this variable likely leading to more runoff and cDOM input) and zooplankton biomass (Fig. 4A). Increasing zooplankton biomass was related to negative effects on Glenodinium and Dinobryon and a positive effect on Trachelomonas (Fig. 4A). However, when only lakes with extreme higher or lower zooplankton biomass values were examined through NMDS-PERMANOVA; no significant differences in overall composition of mixotrophs were found (Supporting Information Fig. S4). This discrepancy could potentially indicate that the effects are nuanced species effects that may be detectable only across the entire zooplankton biomass gradient. Possibly, more global predation effects could have been disentangled with greater temporal resolution (allowing for predator-prey time lags) or after accounting for other unmeasured environmental variables. However, the observation that zooplankton seem to prefer certain types of mixotrophs over autotrophs is in line with findings of Katechakis et al. (2005) that certain mixotrophs may represent a superior food source while others may induce toxic effects on zooplankton.
The inferred preference for mixotrophs by zooplankton suggests that the mixotrophic pathway may help sustain a constant trophic nutrient flow throughout the food web, even as autotrophs become light limited. This would suggest that mixotrophs play a vital role in sustaining lake ecosystems across changing conditions, even as a lake goes from net autotrophic to net heterotrophic. As such, our study supports that a significant fraction of the zooplankton community benefits or at least would not be penalized from allochthonous DOM inputs into lakes (Salonen and Hammar 1986;Bowszys et al. 2014;Berggren et al. 2015;Paczkowska et al. 2017) due to a constant presence of mixotrophy. In fact, Berggren et al. (2014) by using a subset of the same lakes as in the present study found significant evidence for allochthony in the zooplankton communities using multiple stable isotopes. Our analyses reinforce the idea that a significant part of the allochthonous food web pathway likely includes mixotrophs as an intermediate and quantitatively significant trophic link.

Limitations
A limitation of this study is the snapshot sampling approach used. This means that critical factors influencing our response variables may no longer have been present in the system at the time of sampling. Where time lags and a temporal decoupling between different variables are important, determining directional cause-effect relationships is difficult, and can only be done by engaging other in-depth knowledge of lake ecosystem mechanisms, structures, and functions. Although this is a weakness of this dataset, we do not attempt to predict influences of long-term trends in environmental factors. As such, the number of samples, and spatial scale in combination with the extended periods of time between first and last sample, means that it is possible to treat them as independent observations in the context of this study to identify informative relationships.
In our study, as mentioned, we did not measure rates of mixotrophy directly. Instead, phytoplankton taxa with a high capacity for mixotrophy were determined from a literature review (Longhi and Beisner 2010). Furthermore, while recent evidence suggests that some chlorophytes may be capable of bacterivory (e.g., Anderson et al. (2018)), this has not generally been established for lake communities (e.g., Medina-Sánchez et al. (2005)) and we thus considered this group as autotrophic. Ideally, real-time estimates of mixotrophy would be ideal to further this work and to further identify environmental gradients and factors associated with such activity in phytoplankton communities.
Another caveat of this study is the lack of data on heterotrophic microorganisms. Unfortunately, data on heterotrophic activity were not gathered in the initial sampling program and cannot be estimated in any other way. This means that we have no way of inferring how mixotrophs fare in direct competition with heterotrophs under varying conditions. Such data should be examined in future studies of this nature and should shed insight into the responses of mixotrophs relative to heterotrophs especially under net heterotrophic lake conditions.

Conclusions
Success of mixotrophs, as defined by in lake proportion of mixotrophy, is related to the fDOM_4 fraction leading to increased light limitation of the autotrophic phytoplankton community in boreal lakes. Furthermore, in boreal lakes, dissolved CO 2 can be used as a proxy for lake autotrophy to heterotrophy status, whereby a constant mixotroph community maintains an almost constant biomass in both autotrophic and heterotrophic lakes. However, mixotroph biomass appears to become limited when both light and prey abundance are low, at intermediate CO 2 concentrations. Finally, mixotrophic success in our boreal lakes is also driven by a negative impact of zooplankton predation, likely owing to selection of mixotrophic prey which should provide a more stable source of nutrients over autotrophic phytoplankton, with more variable stoichiometry. Thus, for the first time, we demonstrate evidence that mixotrophs are likely a preferred food source for zooplankton grazers under natural conditions at least in a large number of North American boreal lakes.