EARLY SPRING PHYTOPLANKTON DYNAMICS IN THE SUBPOLAR NORTH ATLANTIC: THE INFLUENCE OF HETEROTROPHIC-PROTIST HERBIVORY

To assess the importance of herbivory by heterotrophic protists in relation to mixed-layer depth prior to the spring phytoplankton bloom, we measured phytoplankton growth and heterotrophic-protist grazing rates during the March/April 2012 EuroBasin Deep Convection cruise in the subpolar North Atlantic. We performed 15 dilution experiments during 2-4 visits at one shelf (160 m) and two deep (~1300 m) stations. Of the two deep stations, one had a mean mixed-layer depth of 476 m, whereas the other was stratified (46 m). Euphotic depth averaged ~70 m at both stations. Initial chlorophyll-a varied from 0.2 to 1.9 μg L at the deep mixed layer station and from 0.5 to 1.0 μg L at the stratified station. In 80 % of the experiments, growth rates exceeded grazing mortality rates, regardless of mixed layer depth. Large mixed layer depth coincided with phytoplankton growth and grazing mortality rates that varied over a similar range from ≤0 to 0.6 d, and to an average grazing-impact representing 50% of primary production (PP). At the stratified station, phytoplankton growth rates varied from 0.18 to 0.41 d, grazing mortality rates varied from 0.11 to 0.34 d, and a temporal shift from a positive to a negative balance between growth and grazing rates caused the proportion of PP consumed to increase from 60% to 180%. Variations in in situ chlorophyll-a could not be explained where the mixed layer was deep, whereas at the stratified station the balance between rate estimates of phytoplankton growth and grazing mortality rates explained 98 % of measured changes in chlorophyll-a. These results suggest a difference in the dominant surface loss process at the two stations: grazing at the stratified station vs. potential sinking aided by vertical mixing where mixed layer was deep.


Introduction
In the subpolar North Atlantic, the yearly cycle of primary production (PP) is dominated by the annual recurrence of the spring phytoplankton bloom. The seasonal increase in phytoplankton biomass is of large ecological significance: not only does the bloom fuel marine food webs, it also influences earth's climate, since the associated drawdown of atmospheric CO 2 and the consequent sinking of some of the bloom biomass (Buesseler et al. 2007, Sarthou 2005, Turner 2002) contribute substantially to the strength of the global biological pump (Takahashi et al. 2009).
For a bloom (i.e. an accumulation of biomass) to occur, net phytoplankton population growth rate (i.e. accumulation rate) needs to be positive, that is phytoplankton intrinsic growth rate has to exceed the rate at which production is lost (Banse 1994). The process can be described by the equation r = µ -l, where r is the phytoplankton biomass accumulation rate, µ is the phytoplankton growth rate, and l represents the rate of phytoplankton losses (Behrenfeld 2010).
From the earliest days of the extensive research devoted to identifying what triggers the North Atlantic spring bloom, of the two terms involved in the equation describing a bloom, µ has received the most attention (Behrenfeld & Boss 2014). In particular a large focus has been placed on the influence on µ of one physical variable: mixed layer depth (MLD), a proxy for, yet not always representative of, the actively mixing layer (Ferrari et al. 2014, Taylor & Ferrari 2011. Starting with Sverdrup (1953), the idea that the North Atlantic spring bloom begins when the mixed layer shoals above a "critical depth", i.e. the depth of a mixed layer within which integrated phytoplankton production and losses are equal, has served as a paradigm in the understanding of bloom formation (Henson et al. 2006, Siegel et al. 2002, Sverdrup 1953. Numerous observations, however, have been reported of early spring surface increases in phytoplankton biomass preceding stratification (e.g. Dale et al. 1999, Townsend et al. 1992,1994, challenging Sverdrup's classical bloom model. Yet consequent new hypotheses have continued to focus on potential factors driving µ, all involving the extent of vertical mixing: for example, rates of turbulent mixing (Huisman et al. 1999(Huisman et al. , 2002, heat-flux induced weakening of turbulent mixing (Ferrari et al. 2014, Taylor & Ferrari 2011, and eddy-driven stratification (Mahadevan et al. 2012). Thus traditionally, a disproportionate emphasis has been given to µ with the loss term being less studied.
Of all losses affecting PP, the largest is due to grazing (Banse 1994). In particular, herbivory by ubiquitous <200 µm heterotrophic protists (HP), such as ciliates and dinoflagellates, has been identified as the major fate of ocean PP (Calbet & Landry 2004, Strom 2002). Thanks to their diverse feeding strategies, protist grazers can access a broad range of prey sizes, from bacteria to prey larger than they are (Aberle et al. 2007, Sherr & Sherr 2002. HP grow at rates similar to the cells they eat, allowing their numbers to often increase quickly after an increase in available prey (Sherr et al. 2003). From a plethora of studies performed across oceans to measure HP grazing rates, HP grazing impact has been estimated to average ~69% of PP (Calbet & Landry 2004, Schmocker et al. 2013. Although temporal and spatial exceptions exist, in which other loss processes such as viral lysis (Brussaard 2004) or nutrient starvation ) control phytoplankton biomass, HP herbivory has been established as the most significant loss factor in PP.
From the research that has considered the role of grazing losses in phytoplankton blooms, a consensus has emerged that seasonal high-latitude blooms happen because grazing cannot keep pace with phytoplankton growth. Various mechanisms have been considered, including proposed phytoplankton predationavoidance strategies (Irigoien et al. 2005), and constraints on HP growth rates at temperatures <5° C (Rose & Caron 2007), or low prey pre-bloom availability , Sherr et al. 2013, all yielding µ in excess of grazing.
Recent work by Behrenfeld and colleagues (Behrenfeld 2010, Behrenfeld & Boss 2014, Behrenfeld et al. 2013) has re-examined the importance of the physics of MLD, by considering its effects not only on phytoplankton growth as has been traditionally done, but also on the magnitude of grazing pressure. Behrenfeld (2010) suggested that a key process influencing variations in the North Atlantic phytoplankton biomass is the alteration by vertical mixing of the balance between µ and grazing. According to Behrenfeld's "dilution-recoupling" hypothesis (Behrenfeld 2010), deepening of the mixed layer in winter reduces predator-prey encounters, decreasing phytoplankton grazing losses below the very low but positive rates of winter phytoplankton growth, thus allowing blooms to initiate during winter. Behrenfeld (2010) further postulated that the gradual seasonal shoaling of the mixed layer "re-couples" predators with their prey, resulting in increased grazing pressure, which curbs phytoplankton biomass accumulation rate.
Despite the potential importance of HP grazing, its role in bloom development remains theoretical, as our understanding of pre-bloom grazing dynamics suffers from a shortage of available in situ grazing rates measured before or during early bloom development. In particular, for the open North Atlantic ocean at high latitudes above 50° N, where winter mixed layer is typically large due to convection (Backhaus et al. 2003), existing bloom-related in situ measurements of HP grazing rates come from studies conducted during or after the bloom , Gaul & Antia 2001, Gifford et al. 1995, Stelfox-Widdicombe et al. 2000, Wolfe et al. 2000, and to our best knowledge, there are no empirical data of grazing rates for the critical period that precedes the bloom.
The present research was performed during the early spring 2012 EuroBasin program "Deep Convection" research cruise, which intended to evaluate the response of the subpolar North Atlantic ecosystem to physical forcing during the transitional period when winter convection gradually weakens (Backhaus et al. 2003 We found that at the oceanic sites, rates of both phytoplankton growth and HP grazing could be substantial, even when the mixed layer was deep, yet in most experiments growth rates exceeded grazing rates, regardless of MLD. Further analysis suggested that grazing was a dominant factor controlling phytoplankton biomass only when the mixed layer was shallow, and not when it was deep.

Sampling sites
Heterotrophic-protist herbivory was quantified during two to four visits at two ~1300-m deep open-ocean sites located in the Iceland Basin and the Norwegian Sea (S1 & S2), and at a ~160-m deep site located on the Shetland shelf (S3) (Fig. 1).

Dilution experiments
We measured HP grazing rates in 15 separate experiments using the Landry & Hasset (1982)  was affixed in order to screen out larger grazers ). We further refer to this < 200-µm fraction as whole seawater (WSW).
For nine experiments (one per visit at each station), we prepared five target dilutions (9, 18, 37, 75, 100 % WSW), which were distributed so as to increase our ability to detect potential non-linearity related to feeding thresholds at low prey concentration, their specific values arbitrarily chosen to facilitate measuring the volumes to be combined. For the first and second experiments at S1, the actual number of dilution levels achieved was reduced from five to four, due to loss of the Bottles were incubated for 24 hours. All incubations took place in on-deck ~250-L tanks. Bottles were suspended mid-water by strapping them onto bungee cords loosely stretched across the length of the tanks, which together with ship motion provided gentle agitation. Incubations were maintained at in situ temperature by flowthrough of ambient seawater. Incubation temperature was recorded at 30-mn intervals using in-tank Hobo data loggers. Incubation temperature was on average 0.9 (± 1.1) °C higher than the temperature at collection depth, however departure occurred mainly during the first leg of the cruise, when differences were the largest at the 1 st experiments at S2 and S3.
To minimize chlorophyll bleaching, which is known to occur in light-sensitive polar phytoplankton (Caron et al. 2000, Smith & Sakshaug 1990, bottles were incubated in black neutral-density mesh-bags that reduced the light to 30% of surface irradiance. Incubations carried at collection-depth irradiance fail to truly replicate the average light regime experienced by cells in a mixed layer (Ross et al. 2011), therefore, in general, the same mesh screen was used regardless of water collectiondepth. However, to investigate the effect of light (31 March) and of collection depth (14 April) on rate magnitudes, a set of two experiments were incubated simultaneously, one with and one without a mesh-bag (See Appendix).

Phytoplankton growth and grazing mortality estimates
Phytoplankton growth and HP grazing mortality rates were estimated from changes in extracted chlorophyll-a (chl-a) over the incubation period (Landry &Hasset 1982). Initial and final chl-a concentrations were determined from triplicate subsamples of each dilution stock and of each replicate bottle respectively.
Subsamples ranged in volume from 60 to 500 ml depending on the in situ chl-a concentration and the dilution level. Chl-a extraction and determination followed Graff & Rynearson (2011), except that extraction took place at room temperature for 12-15 hours in 96% ethanol (Jespersen & Christoffersen 1987).
Apparent phytoplankton growth (k, d -1 ) in each bottle was estimated using the equation k = 1/t ln (P t -P 0 ), where t = incubation time in days, and P t and P 0 are respectively the final and the initial chl-a concentrations. We used these estimates of k to determine the instantaneous phytoplankton growth rate (µ, d -1 ) and the instantaneous grazing rate (g, d -1 ) using two methods.
First, as is customary to the dilution method (Landry & Hasset 1982), the rates were determined from the linear regression analysis of k as a function of the dilution factor, where µ is the y-intercept and g the negative slope of the line. We tested the hypothesis that g = 0 for each regression. We applied dilution factors as determined from measured initial chl-a concentrations in the dilutions, which was on average 1.8 % lower than the target (± 3.9 %). For the first experiment at S3, malfunctioning of the fluorometer yielded inaccurate measurements of initial chl-a concentration in the diluted treatments, thus initial chl-a was assumed to equal WSW chl-a multiplied by the target dilution factor. We found no significant difference between k in undiluted treatments with and without nutrients (two-tailed paired t-test; p = 0.63, 0.21, and 0.15 for station 1, 2, and 3 respectively), and consequently combined all undiluted replicates in our analysis.
One of the dilution method's major assumptions is that k is linearly related to the dilution factor. In case of deviations from linearity, a linear regression is thus an inadequate method to determine g and µ (Worden & Binder 2003). We therefore tested whether the linearity assumption held for all dilution experiments or if significant deviations existed (ANOVA, α = 0.05; Zar 2010). We found significant deviations from linearity in three regular ( we did not find the two-point estimated rates to significantly differ from rates obtained from regression analysis (two-tailed paired t-test, p= 0.10 for µ, p= 0.84 for g). To insure internal consistency, the two-point method was used for all experiments, and rates reported herein are those thus derived ( Table 2).
The grazing impact of HP in terms of the proportion of primary production (PP) consumed was calculated as % PP = g : µ × 100 following Calbet & Landry (2004). For all calculations, negative growth rate and negative grazing rate estimates were corrected to +0.01 d -1 and zero respectively. No % PP was calculated for experiments in which no significant phytoplankton growth was measured. For each experiment, we also calculated biomass-specific grazing rates on phytoplankton (G HP ) using the equation G HP = [(g)(Chl)(C:chl)] / HP, in which g is the estimated specific HP grazing rate, Chl is experiment initial chl-a concentration, and HP is HP biomass (Strom et al. 2007, Strom & Frederickson 2008. As we did not measure C:chl ratio during this study, we used a ratio of 21, a value found to be a good estimate for phytoplankton communities of Norwegian coastal waters (Bratbak et al. 2011), although we recognize that this ratio is highly variable and poorly constrained (Geider 1987, Sathyendranath et al. 2009).
As a mean to assess the importance of the grazing loss term in determining phytoplankton biomass dynamics, using chl-a as a proxy for biomass, we compared observed (i.e. in situ) chl-a accumulation rates (r obs ) to the accumulation rates inferred from the balance between experimentally determined phytoplankton growth and mortality rates (r calc = µ -g). The observed accumulation rate was determined using the equation r obs = 1/t * ln (P t -P 0 ) where P t and P 0 are chl-a concentrations (from initial experiment samples) at the end and the beginning of the time interval t separating two consecutive experiments at the same station. Diatoms, dinoflagellates, and ciliates were enumerated by settling 50 ml for a minimum of 24 h following the Ütermohl (1958) method. Diatoms were identified to genus following Throndsen et al. (2007) and Kraberg et al. (2010). Dinoflagellates were divided into thecate and athecate groups, and when possible further identified to genus following Dodge (1982), or assigned to a morphotype (based on similarity of shape). Preservation of samples in acid Lugol's fixative does not allow differentiating between auto-and heterotrophic dinoflagellates. Furthermore, many autotrophic dinoflagellates are also phagotrophic (Stoecker 1999), thus all dinoflagellates were assumed heterotrophic. Ciliates were divided into loricate (tintinnids) and aloricate groups. Higher taxonomic identification of aloricate ciliates relying on shape was attempted following Strüder- Kypke et al. (2002), to provide a qualitative description of the ciliate community, but due to its lack of reliability (Montagnes & Lynn 1991), it was not used for quantitative analysis. An exception was however made for the obligate mixotroph species Mesodinium rubrum (= Myrionecta rubra; Hansen et al. 2012), which could easily be distinguished from and were not included with other aloricate ciliates.

Biomass and composition of the plankton community
Linear dimensions were measured using ImageJ software (National Institute of Health) from images taken of all dinoflagellates and ciliates contained in each sample and, depending on abundance, of all or a subset of diatom cells (30-300 cells per genus). Cell volumes were calculated from linear dimensions using appropriate geometric shape algorithms.

Ancillary data
We used hydrological data collected by CTD to characterize in situ conditions at the depth of sample collection, as well as general environmental conditions encountered during the study. A total of 20, 14, and 9 full-depth CTD casts were available for S1, S2, and S3 respectively, which we used to generate estimates of mixed-layer depth (MLD), average mixed-layer temperature (T) and salinity, and MLD integrated chl-a concentration. In estimating MLD, we adopted a T threshold criterion of -0.2 o C from a reference depth of 10 m (de Boyer Montégut et al. 2004).
We also used CTD data to estimate surface photosynthetically available radiation (PAR), and depth of the euphotic zone (Z eu ). Due to lack of PAR data from <5 m depth, we estimated surface PAR as the y-intercept (depth = 0 m) of a linear regression of the natural log of PAR profiles, the slope of which yields the coefficient of vertical light extinction (k i ). We then used these estimates of k i to determine Z eu .
We adopted the commonly used definition of Z eu as the depth receiving 1% of surface irradiance (Margalef 1978, Reynolds 2006).

Statistical analyses
Principal component analysis (PCA) was used to summarize environmental variability. Included in the analysis were data from CTD casts used to collect water for the experiments of in situ temperature and salinity, and estimates of MLD, and Z eu .
Before analysis, non-normally distributed data were log-transformed, and to place all variables on a comparable dimensionless scale, data were normalized to a mean of zero and a standard deviation of 1.
Patterns in the composition of the diatom and of the HP assemblages were investigated using the non-parametric multivariate statistics package Primer-E Plankton biomass-based similarity matrices were also used to further compare plankton assemblages based on location (3 levels: S1, S2, S3) and on grazing magnitude. The latter factor was partitioned into three different levels of grazing activity relative to the overall average (zero, below average, and above average). To examine if species composition influenced whether grazing occurred at all, the analysis was repeated using only two levels (0= no grazing, 1= grazing).
To assess the nature and strength of relationships between species composition samples, analysis of similarities (ANOSIM) was performed on biomass-based resemblance matrices. ANOSIM is a non-parametric permutation procedure that computes the global R statistic, which can range from -1 to 1, although negative values are unlikely (Clarke & Warwick 2001). Values approaching 1 indicate greater similarities within a group than among groups, whereas values approaching zero indicate equal similarities within and among groups (e.g. no group associations/clustering).
Additionally, a series of univariate analyses (linear regression and Pearson correlation) were performed using SigmaPlot ® software to examine relationships between grazing rates and a series of potential driving factors. All statistical analyses were performed at an alpha level of 0.05. All rates and other estimates are expressed ± one standard deviation from the mean.

Spatial patterns of in situ conditions
Contrasting environmental conditions over the sampling period distinguished the three stations. Data from all CTD casts performed over the entire duration of the cruise provided evidence that temporal variation in physical parameters were greater among than within stations.
Stations significantly differed in MLD (ANOVA, p < 0.001). S1 had a deep mixed layer, which averaged 476 ± 149 m and showed no shoaling progression. MLD at S2 was one order of magnitude shallower than at S1, averaging 46 ± 16 m. At the shallow (160 m) shelf station S3, MLD always reached the bottom. Consequently, MLD at S1 was repeatedly deeper and MLD at S2 was generally shallower than the euphotic depth, which was estimated to average 70 m (± 18 and 10 m at S1 and S2 respectively). On the shelf at S3, MLD was always deeper than the mean euphotic depth of 50 ± 10 m.
Stations also differed in temperature (T) and salinity. Mixed-layer average T was warmest at S1, where over the sampling period, it averaged 8.61 (±0.23) °C, and Accordingly, conditions recorded in situ at the time of sampling (Table 1) were characterized by significant differences among stations (ANOSIM global R= 0.796, p = 0.002) driven by differences in in situ T and MLD (Fig. 2). Together the first two axes of the PCA explained 89.6% of the variance of the in situ data.

Species composition of the plankton assemblage
There were clear spatial differences in the species composition of the plankton assemblage of S1 and S2, as station-specific samples taken over a period of 33 days resembled each other most (Fig. 3). Both the diatom (Fig. 3a) and the HP assemblages ( Fig. 3b) from each sample were strongly associated with location (ANOSIM p ≤ 0.002), and S1 and S2 differed the most (p = 0.002). Temporal variability of HP assemblage among station-specific samples was greater at S1 than at S2, whereas the reverse was true for diatoms, which at S2 were scarce (see below). At S3, a shift in phytoplankton species composition and biomass (see below) caused the diatom assemblage in the two experiments to be <40% similar. Corresponding HP assemblage samples were <50% similar, with the sample from the 1 st experiment at S3 resembling those of S2 the most. Both the diatom and the HP assemblages correlated with the multivariate pattern of environmental data characterized by the PCA (RELATE Spearman correlation= 0.518 and 0.47 respectively, p= 0.002).
Based on a C:chl-a ratio of 21 (Bratbak et al. 2011), diatoms would have comprised 9-40% of total chl-a, increasing to their maximum on April 10 and decreasing thereafter.
The genera Chaetoceros and Pseudonitzschia dominated the diatom community, respectively representing up to 86 % (April 9) and 62% (April 28) of total diatom biomass. In contrast at S2, diatoms were quasi-absent (average biomass <0.1 µg C L -1 ), and the >10-µm chl-a fraction never exceeded 8 % of total chl-a (M. Paulsen pers. comm.). Although not included in estimates of autotrophic biomass, a large number of cryptophytes identified as Teleaulax spp. were present in the S2 samples. At S3, diatoms increased in biomass from 0.56 to 87.7 µg C L -1 between the 1 st and the 2 nd visits (Table 3). At the 2 nd visit, the species Detonula pumila and Ditylum brightwellii together constituted most of the diatom biomass (64 % and 25 % respectively).
There was no within-station correlation between heterotrophic biomass and chl-a concentration (Pearson correlation, S1 & S2 p= 0.83). One concern was that collection depth, which differed among experiments, might have affected concentration of protistan grazers and by extension grazing rates, but it did not significantly influence either their numerical abundance (p>0.45) or their biomass (p>0.43).

Results of dilution experiments
Over the entire sampling period and across all stations, phytoplankton growth rates ranged from -0.06 to 0.63 d -1 and mortality rates due to HP grazing ranged from 0 to 0.56 d -1 (Table 2). In all but three measurements, growth rates exceeded grazing mortality rates (Fig. 5). The magnitude and variability of growth and grazing rates at S1 and S2 differed, with S1 exhibiting both higher rates and higher variability. At S1 growth and grazing rates varied over the same range (0 to 0.6 d -1 ), although average growth rate 0.35 (± 0.03) exceeded average grazing rate 0.25 (± 0.04) d -1 (Table 2).
There was one exception to the general decoupling between growth and grazing rates at S1: on April 10, rates were highly coupled (0.60 and 0.56 d -1 respectively), and corresponded to the highest initial concentration of chl-a (1.9 µg L -1 ) of all experiments (Table 2).
At S2, growth rates ranged from 0.18 to 0.41 d -1 and grazing rates ranged from 0.11 to 0.34 d -1 . Growth and grazing rates had similar averages (0.24 ± 0.02 d -1 and 0.22 ± 0.03 d -1 respectively) ( Table 2). On the last two sampling dates, the balance between phytoplankton growth and grazing rates changed from positive to negative.
On the Norwegian shelf (S3), only two experiments were performed at a twoweek interval. The first experiment yielded no detectable grazing, and a very low grazing rate (0.04 d -1 ) was measured the second time (Table 2), whereas phytoplankton growth rates were similar on both dates (0.23 and 0.27 d -1 ).

Influence of grazing on dynamics of phytoplankton biomass
The two oceanic stations differed in the level to which in situ chl-a variations followed the dynamics inferred from the rates. Based on the average balance between growth and grazing rates and assuming no other losses than grazing, phytoplankton population at S1 would (on average) have doubled approx. every week, whereas at S2, it would have doubled approx. every month. At S1, measured variations in chl-a did not match those inferred by the rate estimates (R 2 = 0.10, p= 0.61) (Fig. 7a). We measured a 10-fold increase from 0.2 to 1.9 µg L -1 between March 26 and April 10, which clearly exceeded the ~zero growth rates measured in the first two experiments.
Initial experimental chl-a remained ~1 µg L -1 for the rest of the sampling period despite growth rates exceeding grazing rates. CTD profiles, however, show various sub-surface chl-a increases at S1 beyond the maxima measured in our experiments (e.g. Fig. 8), which corresponded to small temperature increases, and after which (from April18) the vertical extent of chl-a was well below the euphotic zone and closely coincided with MLD (Fig. 8). Based on CTD data, MLD-integrated chl-a concentration increased from ~40 mg m -2 at the 1 st visit to 230-250 mg m -2 during visits 3 and 4.
In contrast with S1, at S2 the observed variation in chl-a closely matched the balance between experimentally estimated rates (R 2 = 0.98, p= 0.009) (Fig. 7b).
During the 1 st half of April, phytoplankton growth rates exceeded grazing rates and surface chl-a doubled from a mean of ~0.5 to ~1.0 µg L -1 . For the last two experiments at the end of April, grazing rates exceeded phytoplankton growth rates.
Coinciding with the decoupling of growth and grazing rates on April 23 (balance = -0.2 d -1 ), an overnight 20 % decrease in chl-a was observed. This decrease also corresponded to an overnight change in MLD from 29 m to 68 m.
At S3, based on estimated rates and assuming no other losses than grazing, phytoplankton population would have doubled every ~3 days, twice more often than indicated by the two-week increase in chl-a from 0.5 to 2.7 µg L -1 . This 5-fold increase coincided with a tripling of HP biomass. The increase in HP biomass corresponded to a 31 % decrease in grazer abundance and thus to a shift to a four times larger average size of grazers.

Few specific drivers of protistan herbivory
Several factors potentially governing the magnitude of grazing rates were examined. Experiments that yielded no detectable grazing were not included in univariate analyses. No correlation existed between grazing magnitude and either in situ or incubation temperature (p≥ 0.85). One of the parameter of interest in this study was MLD, but we found no within-station significant correlation between MLD and grazing rates (Pearson, p≥ 0.23). At S2 however, the lowest grazing rate (0.11 d -1 ) was measured on April 13, when unlike all other times, the sample collection depth (35 m) was greater than our estimate of MLD (30 m), and when ciliate biomass was the lowest (2.3 µg C L -1 ). Furthermore, maximum S2 grazing rate (0.34 d -1 ) was measured when MLD was the shallowest (29 m). For S1 and S2 combined, grazing rates decreased with increasing collection depth (R 2 = 0.54, p= 0.016), but within-station relationship was not significant (R 2 = 0.72 and 0.34 and p = 0.069 and 0.304 for S1 and S2 respectively). When combined, S1 and S2 grazing rates significantly correlated with chl-a (Pearson coefficient= 0.721, p= 0.019), but no significant correlation existed within each station (p≥ 0.18). No significant correlation existed between grazing rates and either HP biomass (p= 0.94 and 0.08 for S1 and S2 respectively), or HP numerical abundance (p> 0.58). At S2, daily grazing rates transformed into carbon ingested per day using a C:Chl ratio of 21 tended to covary with HP biomass (Fig. 6) but the correlation was not statistically significant (rho = 0.837, p = 0.077). We further investigated the effect of species composition of each the autotrophic and the heterotrophic assemblage on grazing patterns/rates, using indices of grazing rate magnitude as a factor in analyses (see methods) and found no significant correlation.

Discussion
Our study is, to our best knowledge, the first among a plethora of published field measurements (e.g. Calbet & Landry 2004, Schmoker et al. 2013 to provide estimates of pre-bloom heterotrophic protist grazing rates in the subpolar North Atlantic. Such rare estimates are much needed to complement proposed hypotheses , Rose & Caron 2007, Irigoien et al. 2005) about the role of HP feeding in the development of phytoplankton blooms. We also examined how MLD may have modulated the balance between µ and g, a process that has been proposed to be a major factor controlling variations in phytoplankton biomass, including when the spring bloom initiates (Behrenfeld 2010).
In this study, significant, positive rates of protistan herbivory were often measured, representing a potentially substantial loss of phytoplankton biomass, yet for the most part grazing could not keep pace with phytoplankton growth, regardless of MLD. Our findings further suggest that, at the two open ocean sites with contrasting MLD, different processes were driving phytoplankton losses from the surface layers.

HP grazing rates and grazing impact
Major assumptions of the dilution method and deviations thereof have been discussed at length (Dolan et al. 2000, Moigis 2006, Agis et al. 2007, and others summarized in Schmoker et al. 2013) and are not addressed here. More rarely mentioned is that the dilution method assumes that mortality rates are entirely due to grazing, when they may include phytoplankton mortality due to physiological senescence (Franklin et al. 2006), and perhaps more importantly viral lysis. Although the magnitude of virus-induced mortality is poorly constrained due to the limitations of available methods, it may be significant (Brussaard 2004), potentially varying with the type of trophic interactions regulating carbon flow within the plankton (Ory et al. 2010). Due to their influence on carbon and energy fluxes (Moore et al. 2004, Suttle 2007), these two processes would deserve to be more routinely evaluated.
This being acknowledged, our results indicate that HP collected from F-max were active grazers of chl-a at the two oceanic sites sampled during this study, consuming 26-94 % and 45-242 % of daily PP at S1 and S2 respectively, i.e. >60% in 9 out of 11 experiments that yielded >0 grazing. With few exceptions, these values are similar or greater than the average estimate for other oceanic regions (70%) or polar and temperate regions (60%) (Calbet & Landry 2004). Grazing impact was however variable, as is characteristic of most studies, including previous studies conducted in the region at different times spanning May to July (Table 4).
In contrast to the % PP consumed, our estimates of grazing rates, which ranged from 0 to 0.56 d -1 , were at the lower end of the range of rates (0-1.48 d -1 ) measured in previous studies (Table 4). Average grazing rates at the two oceanic stations differed slightly (by 0.03 d -1 ) with a maximum difference of 0.2 d -1 . When comparing the two sites, temperature (T) difference has to be considered, since T influences ingestion rates (Hansen et al. 1997, Verity et al. 2002. Based on published values of Q 10 for ingestion rates (Verity et al. 2002), the largest difference in T between the two oceanic stations (2°C) could have produced a ~0.1 d -1 difference in grazing rates, thus T may have marginally impacted the magnitude of the rates measured.
Grazing rates were particularly variable at S1. HP distribution has been found to be patchy at very fine scales (Montagnes et al. 1999), which may inherently confer variability among grazing rates resulting from different experiments due to differences in the species and size composition of HP assemblages. The relationship between variability of HP assemblage and grazing rates seem to be corroborated by our observations that HP assemblages were more variable at S1 than at S2. were unusually coupled, results may have been affected by horizontal advection of a plankton patch. Treating our results from April 10 as an "outlier" would bring a temporal pattern to the rates we measured, with the first phase of the sampling period being associated with ≤ 0.1 d -1 growth and grazing rates, followed by a steady increase to their highest values on the last sampling date. We were not, however, able to firmly establish if sampling occurred in a different water mass.

MLD and mechanisms of uncoupling between µ and g
In this study, we wanted to address the question of how much MLD may modulate the balance between µ and g. The "dilution-recoupling" hypothesis proposes that winter mixed layer deepening dilutes predators and prey, reducing grazing and causing accumulation rates to become positive, whereas "recoupling" occurs once the mixed layer deepening stops, and eventually biomass accumulation rates become negative when vernal stratification sets in (Behrenfeld 2010). Most of the experiments we conducted yielded growth rates that exceeded grazing rates, providing a potential mechanism for phytoplankton biomass to accumulate and potentially form a bloom, regardless of mixed layer depth. MLD did seem to influence the magnitude of the balance between µ and g, which was larger at S1 than at S2. Nevertheless, although the magnitude of the balance is important in setting the accumulation rate, it is its sign (>0 or <0) that ultimately controls the potential for biomass to accumulate. The majority of our results indicate that MLD was not a main determinant of whether the µ-g balance was positive or negative. In particular, although there were exceptions, the fact that growth rates exceeded rates of grazing losses, including at S2, where the mixed layer depth was approx. half as deep than the euphotic depth, suggests that stratification by itself may not always be sufficient for grazing to become large enough to decrease phytoplankton biomass. At S2, early stratification had occurred.
Although its mode of formation has not been firmly established, it was unlikely related to vernal warming, as air temperature ranged between 2-5 °C. Despite of stratification, food and/or temperature conditions at the end of winter may not have been conducive to active growth of HP, at least in the first part of the sampling season, when growth rates exceeded grazing rates.  presented evidence that HP cannot prevent blooms due to generally low food availability, leading to low growth rates and low biomass of HP under non-bloom conditions. We may have underestimated HP biomass by not including heterotrophic nanoflagellates (HNAN), which can contribute significantly to HP biomass (Stelfox-Widdicombe et al. 2000). According to our estimates, HP biomass remained low for the entire sampling period at all stations, particularly at S1 where it corresponded to the lower limit of the range of values previously published for the region at other times of the year (Table 4). Our statistical analyses did not provide evidence of a direct relationship between HP biomass and grazing rates. To further examine this relationship, we calculated biomass-specific grazing rates (G HP ) and compared them to maximum laboratory-determined rates. Published estimates of G HP generated by laboratory experiments reported in Hansen et al. (1997)  only once. At all other times, G HP values were < 3 d -1 at S1 and < 1 d -1 S2, and even when adjusted for temperature, these values remain low, especially at S2. Overall low values of G HP could have resulted from the C:Chl ratio of 21 used in their estimates.
Although low, in our study this ratio may have been representative of the physiological state of phytoplankton having to acclimate to early spring low levels of light and active vertical mixing in abundant nutrients, all of which would contribute to substantial cell -1 chl-a (Geider1987). Using this ratio to determine the relative contribution of diatoms to total chl-a yielded estimates that adequately compared to estimates made by other investigators on the cruise (M. Paulsen, pers. comm.).
Furthermore, higher estimates of HP biomass (by including HNAN) would further reduce G HP values. Low values indicate that HP were feeding below their potential rates. Although the smallest grazers may have been feeding exclusively or alternatively on bacteria, generally low G HP suggest HP may have been food limited.
For the majority of the experiments, HP were found to consume between 25 and 100% of their body carbon daily. At S1, however, these values varied, and herbivory represented up to 400% of HP carbon biomass, but values <500% are considered low rations . These values would again indicate that HP were food limited, although feeding on alternative prey such as bacteria or other protists cannot be ruled out.
The uncoupling between µ and losses needed for phytoplankton biomass to accumulate can be achieved if g is kept low or if µ is large/increases. At the time our study took place, mixed layer deepening at S1 had stopped. Instead we had entered the period when hypothesized increases in grazing pressure are compensated by increases in µ in response to improved light conditions (Behrenfeld 2010). Contributing to the uncoupling at S1, were high values of µ: in two thirds of the experiments at S1 growth rates were equivalent to doubling times of 1-2 days. Our estimates of µ were based on changes in chl-a, and thus could have been overestimated if differences in the light regime experienced by the cells in incubation vs. in situ caused the cells to increase their pigment concentration during the incubation period, however such photoacclimation may be too slow to have significantly affected growth rate estimates measured over a 24-hour period (Landry et al. 1995). Increases in surface chl-a recorded in CTD profiles, whether they indicate in situ growth or horizontal advection, do provide evidence of the capacity of phytoplankton to sustain substantial growth rates at the latitudes and time of year we sampled. In fact, such increases in surface phytoplankton biomass prior to stratification are characteristic of the process leading to the bloom climax (Townsend et al. 1994, Backhaus et al. 2003. Furthermore, at S1, a difference in which size fraction of the phytoplankton had highest growth rates and which size was consumed may have maintained the uncoupling between µ and g. In particular, diatoms likely enhanced µ, but may not have been readily consumed. This remains speculative, as we did not measure either µ or g for different size fractions separately. Nonetheless diatoms can grow at high rates (Smayda 1997cited in Tillmann 2004 and are physiologically adapted to the highly variable and low light regime induced by frequent mixing (Weeks et al. 1993). Thus diatoms likely enhanced total phytoplankton growth rates. Interestingly, µ values at S2, where very few diatoms were observed, were in general lower than at S1.
Furthermore, measures of HP grazing on mixed phytoplankton assemblages often show higher grazing rates on small cells (Gifford et al. 1995, Strom et al. 2007). In their investigation of taxon-specific grazing, Gaul & Antia (2001) reported grazing avoidance of diatoms and selective preference for small cells, although in their study, selective grazing may have been driven more by active growth of prey than by prey size. Although HP as a group can graze on a broad range of prey sizes, and individual grazers can adapt their own morphology to the size of the available prey, not all grazers can feed on all sizes (Strom 2002). Heterotrophic dinoflagellates are often the major consumers of diatoms (Lawrence & Menden-Deuer 2012, Levinsen & Nielsen 2002, Sherr et al. 2013, Strom & Frederickson 2008), but they were not abundant in our study. Interestingly, the highest grazing rates measured at S1 were associated with few, but at other time absent, >50 µm dinoflagellates. Dinoflagellates' growth rates, low in comparison with other HP (Tillmann 2004), may have prevented them to keep pace with diatom growth. Ciliates, which typically feed on smaller particles (Strom 2002), may not have been actively consuming diatoms. Although they have been observed to feed on diatoms (Sherr et al. 2013), this type of feeding is likely restricted to larger (>50 µm) types (Aberle et al. 2007), a few of which only appeared at S1 after 4/18, possibly responding to an increase in larger prey. Similarly at S3, dinoflagellate biomass was among the lowest we observed during our study (<0.4 µg L -1 ) and did not increase during the two weeks between sampling dates.
Simultaneously, ciliate size substantially increased concurrently with the change in the phytoplankton size distribution. Yet even the larger ciliates may only have been able to feed on diatoms at a slow rate due to a possible increase in the time needed to handle the prey (Irigoien et al. 2005). Clearly more has to be learned about the relation between the size structure of the phytoplankton community and the prey preferences and feeding interactions of the various predators. The size-related loophole hypothesis proposed by Irigoien et al. (2005) may well apply to early spring phytoplankton dynamics at high latitudes, when diatoms first start to grow.
High variability of the physical environment at S1 may have produced a patchy grazing response. Small but distinct surface increases in T recorded in CTD profiles support the idea that oceanic heat uptake was at times sufficient to confer stability to the water column and stall convective mixing, which even if transient, was sufficient to provide a window of opportunity for growth (Townsend et al. 1994, Taylor & Ferrari 2011. Such periods of increased water column stability were intermittent with periods of deep mixing, as evidenced by the presence of substantial chl-a at large depths. This variability in the convective mixing regime, which is believed to be diurnal (Taylor & Stephens 1993), may have influenced the variability of S1 grazing rates, HP biomass, and HP species composition among experiments.
Thus at S1, a slow response of HP to increases in prey due to the episodic nature of the physical environment, and a seasonal shift in the size structure of the prey field, may have favored the uncoupling between µ and g.

Mixed layer depth and grazing control of phytoplankton biomass
At S1, the balance between µ and g could not account for ambient (i.e. depth of collection) chl-a variations. We are mindful of the caveats and assumptions associated with comparing observed rates of change in chl-a with the balance between experimental estimates of µ and g. Differences between the two may inevitably result from the difficulty to duplicate in incubation experiments all field conditions that can affect chl-a concentrations. This may be particularly true in regions of active deep mixing. In incubation bottles, plankton assemblages are artificially kept at one depth and isolated from the ambient turbulence, which vertically re-distribute plankton cells, thus changing light availability (Ross et al. 2011), and possibly altering encounter rates between predators and prey. Such differences can lead to experimental estimates that vary from true in situ values. Furthermore, our sampling frequency imposed long time intervals between experiments, obviously producing gaps in our data.
While the described caveats may have contributed at S1 to the lack of agreement between ambient changes in chl-a and the balance between µ and g, most previous studies have shown that the balance between phytoplankton growth and HP grazing rates can rarely account for the variability of chl-a, possibly due to a general lack of equilibrium between growth and loss processes (Schmoker 2013). Such observations have been made even when sampling frequency was high (Lawrence & Menden-Deuer 2012). Furthermore, ambient departures from the µ -g balance may be inevitable considering that the dilution method provides estimates of potential grazing rates, which are obtained from a truncated plankton assemblage, from which mesozooplankton grazers, known to feed both on phytoplankton and heterotrophic protists (Calbet & Saiz 2005, Saiz & Calbet 2011, have been removed. The importance of such artifact in dilution experiments remains to be determined (Schmoker et al. 2013).
Nevertheless a poor match between ambient chl-a variability and the balance between our rate estimates would suggest a minimal control of grazing relative to other processes on the dynamics of chl-a. In particular, our data do not support the idea that decreases of surface chl-a, observed at S1 both in our experiments and in chla vertical profiles, were due to grazing, when phytoplankton growth rates exceeded grazing rates. In contrast to surface layers, vertically integrated phytoplankton accumulation rate was overall positive. The presence of chl-a below the euphotic zone, where phytoplankton growth cannot be sustained, indicates that down-mixing was a major loss process of primary production from the surface layers. As Backhaus et al. (2003) justly remarked, if concentrated within a shallow mixed layer, the observed vertically integrated biomass would be similar or even surpass spring bloom concentrations. Additionally, as the mixed layer shoals, some of the phytoplankton will inevitably become trapped below the thermocline (Backaus et al. 2003, Behrenfeld et al. 2013. Thus intermittent deep mixing in early spring may be an important mechanism for export of carbon, which may exceed sinking losses associated with surface blooms, much of which may be respired through grazing (Behrenfeld 2010).
In contrast with S1, at S2 the balance between phytoplankton growth and grazing-mortality rates was a reasonably good predictor of in situ phytoplankton population dynamics. This suggests that losses due to sinking may have been limited, and that the majority of the losses incurred were due to grazing, making grazing an important determinant of variations in phytoplankton biomass. Among phytoplankton species, diatoms are believed to drive carbon export because their heavy silicate frustules cause them to sink (Sarthou et al. 2005, Smayda 1970). At S2 the phytoplankton community was dominated by pico-and nanophytoplankton, and diatoms were rare. Small particles are less likely to sink, and their vertical retention may increase grazing opportunities. Thus the importance of grazing in determining variations in ambient phytoplankton biomass at S2 was likely influenced by species composition of the phytoplankton community.
We could not firmly establish the source of stratification at S2. Mesoscale variabilities are frequent in the Norwegian Sea (Hansen et al. 2010). S2 sat near the Iceland Faroe front, and thus stratification could have resulted from the mix of North Atlantic water transported in the Faroe Current with East Iceland Current water (Hansen & Østerhus 2000), or an eddy could have developed along the front. Early stratification other than through surface warming can also be driven by the formation of eddies induced by the slumping of the North-South density gradient associated with the latitudinal differences in temperature (Mahadevan et al. 2012). Such stratification is believed to trigger early patchy phytoplankton blooms (Mahadevan 2012) dominated by diatoms (Alkire et al. 2012). At S2, however, no bloom and few diatoms were observed during the 30 days of the study, despite availability of ample macronutrients. Low contribution of diatoms to total phytoplankton has previously been observed in a nutrient-rich mesoscale eddy with shallow mixed layer (Stelfox-Widdicombe et al. 2000). Thus eddy-stratification may not always result in PP being dominated by diatoms.

Drivers of grazing magnitude
In this and previous studies (e.g. Lawrence & Menden-Deuer 2012, Strom et al.2007), drivers of grazing magnitude remain elusive. As in previous studies (Lawrence & Menden-Deuer 2012, Menden-Deuer & Frederickson 2010), there was no direct relationship between chl-a and grazing rates, confirming that many factors other than a coarse metric for phytoplankton quantity can affect grazing rates (Sherr & Sherr 2007). Failing to identify specific drivers of grazing magnitude may be due to the fact that the dilution method provides bulk estimates of grazing, which result from a poorly constrained multitude of complex feeding interactions. Planktonic trophic links include mixotrophy and omnivory , Flynn et al. 2012) and trophic cascades (Calbet & Saiz 2013). They involve taxonomically diverse organisms, which span a large size range and exhibit a variety of prey preferences and prey selection , Montagnes et al. 2008, and whose feeding behaviors respond in specific ways to the surrounding physical, chemical, and biological conditions (Caron & Hutchins 2012). Teasing apart planktonic food webs both through laboratory experiments and in the field remains challenging but necessary to increase our ability to predict grazing losses.

Concluding remarks
Our study is one of the first to document HP grazing for the subpolar North Atlantic during the early spring period that precedes the seasonal increase in surface phytoplankton biomass. Although models can resolve the effect of physical and biological forcing on PP much more comprehensively than logistically-intensive field experiments ever will, such models can only be accurately parameterized if field measurements of key rate processes such as grazing are available. Our ability to predict how North Atlantic PP will respond to warming ultimately depends on a better understanding of what controls variations in phytoplankton biomass. Results of several modeling studies (Boyd & Doney 2002, Le Quere et al. 2003, Sarmiento et al. 2004 predict that, at high latitudes, increased thermal stratification will result in greater light and longer growing season afforded to photosynthetic organisms, which should increase present-day light-limited PP (Doney 2006, Riebesell et al. 2009).
Although more studies comparing losses incurred by phytoplankton under different conditions of mixing depths are needed to generalize our findings, the different dominant loss factor -sinking in deep mixed layer, and grazing in shallow mixed layer -suggested by our results may imply that a longer period of stratification could reduce the export of organic carbon that occurs due to deep mixing before stratification and the spring bloom climax, whereas more PP could potentially be lost to respiration associated with HP grazing. In the field, higher geographical and temporal sampling resolution will be needed to capture the dynamics of the biophysical factors driving coupling/decoupling between phytoplankton growth and grazing-mortality rates.

Light and collection depth experiments
This appendix reports and discusses results of two experiments conducted at S2 to investigate the effect of light (3/31) and of collection depth (4/14) on rate magnitudes. These two experiments were conducted in parallel to a regular experiment and similar methods were applied, except that they were incubated at surface irradiance, without the use a mesh-bag.
On March 31, the same bulk assemblage collected from a 20-m depth (night cast) was used in two parallel experiments incubated at two light regimes, in order to assess the effect of light on phytoplankton growth and grazing rates. Phytoplankton growth rate was 35% higher when replicates were incubated at mixed-layer adjusted light (0.34 ± 0.04 d -1 ) than when unprotected by light screening mesh (0.25 ± 0.04 d -1 ).
On April 14, water collected from a 5-m depth was incubated at surface irradiance simultaneously with the experiment using water from F-max on that day.
Both depths had similar in situ chl-a concentration (0.6 µg L -1 ). Heterotrophic biomass was similar in the two samples (7.1 µg C L -1 at 5m and 7.9 µg C L -1 at 30 m).
The decrease in phytoplankton growth rates in response to increased incubation irradiance may underline one of the caveats associated with estimating growth rates using the dilution method based on measured changes of chl-a (Landry et al. 1995). Exposing cells to higher irradiance relative to in situ likely caused phytoplankton cells to photo-acclimate, i.e. adjust their cellular pigment concentration (Geider 1987), yielding artificially lower estimates of growth. We did not estimate total phytoplankton numerical abundance before and after the experiment, so we cannot be certain whether growth rates were lower due to changes in cellular pigment concentration or an actual difference in the doubling time. In the second experiment, in which plankton assemblage were collected at two different depths, differences in growth rates could have resulted from differences in phytoplankton species composition, which we did not verify. Ross et al. (2011) showed, however, that growth rates based on chl-a or on carbon differed considerably from each other, particularly near the surface, due to cells photo-acclimating to higher light intensity by reducing their chl-a synthesis while uptaking carbon and thus increasing their C:Chl ratio. In experiments in which plankton were incubated in mesocosms at two different depths (and thus at two light regimes), chl-a increase at low light was found to be mostly due to photo-acclimation, and phytoplankton growth rates were higher when plankton were exposed to higher light (Calbet et al. 2012). Photo-acclimation, however, should not alter grazing estimates, as long as it affects all dilutions used to compute grazing rate equally (Landry et al. 1995).
In contrast to growth rates, which responded similarly to light in the two experiments, grazing rates responded differently. In the March 31 st experiment, light showed no effect on grazing. Calbet et al. (2012) measured similar grazing rates in differently illuminated mesocosms. On the other hand, Strom (2001) found that light can enhance ingestion, digestion, and growth rates of herbivorous protists, however the enhancement was relative to rates obtained in total darkness. Thus the difference in light between our two treatments may not have been large enough to produce an effect. Further experiments would be necessary to determine a light threshold above which grazing rates significantly increase. The April 14 th experiments contained two variables (depth and incubation light) that confound the interpretation of the results.
Nevertheless, if we take clues from the March 31 st experiment, in which light only affected growth rates and not grazing rates, then the observed difference in grazing rates between the April 14 th treatments were likely due to a difference in collection depth, or to factors associated with it. In the 5-m sample, which yielded the higher grazing rate, heterotrophic protists were twice as abundant but on average half the size of those collected at 30 m, which may have played a role, potentially through an increase in encounter rates. Additionally, since the autotroph community was dominated by small cells, a large proportion of it being composed of picoplankton, grazing rates may have been enhanced because of a better match between the smaller grazers found at 5 m and the size spectrum of the prey.  Table 2. Initial chlorophyll-a concentration (Chl-a, µg L -1 ), phytoplankton growth (µ) and grazing mortality (g) rates, and grazing impact as % of primary production (% PP) consumed (=100 × g/µ). Rates are given per day. Values in parentheses represent standard deviation. In the 2-point analysis (see text for details), rates were determined using the undiluted and the lowest dilution treatments only. Net apparent growth rate in the undiluted treatment (k 1 ) is also given. Significance of deviation from linearity (DL) for linear regression (alpha= 0.05): NS= not significant S= significant. For calculations, <0 phytoplankton growth and grazing rates (marked *) were set to 0.01 and 0 respectively.   Table 4. Results of studies previously conducted at high latitudes of the subpolar North Atlantic during or after the spring bloom to quantify heterotrophic-protist herbivory, including grazing rates (g), proportion of primary production consumed (% PP), chlorophyll-a concentrations, and numerical abundance (10 3 cells L -1 ) and biomass (µg C L -1 ) of heterotrophic protist grazers (HP). Results from the present study are summarized for comparison.