Role of protozooplankton in the diet of North Sea autumn spawning herring (Clupea harengus) larvae

The role of small prey (< 200 µm) in larval marine fish nutrition is largely understudied. Here, we explore the contribution of protozooplankton (PZP 20–200 µm) to larval fish diets, compared to metazoan microzooplankton (MZP 55–200 µm). More specifically, we tested whether the contribution of PZP increased during the low productivity season and decreased as larvae grow. We used North Sea autumn spawning herring (Clupea harengus) as a case study, as it is a key species with high commercial and ecological importance. In autumn and winter, the potential PZP and MZP prey was dominated by cells < 50 µm (mainly Gymnodiniales, Pronoctiluca pelagica,Tripos spp. and Strombidium spp.), while copepod nauplii and copepodites where more abundant in autumn than in winter. Based on their trophic enrichment (∆15 N), larvae preferentially grazed on small MZP < 50 µm rather than PZP both in autumn and winter. Larvae of different body size (range 8–14 mm standard length) fed at the same trophic level but on different prey (similar δ15N but different δ13C). Growth rates (based on RNA/DNA estimates) were similar in autumn and winter, suggesting that growth was not affected by station-specific differences in the composition of the prey field. Our results not only underscore the important role of MZP on larval herring diets both in autumn and wintertime, but also emphasize the limitations of bulk stable isotope analysis. Given the current low recruitment in North Sea herring, these results provide significant information for future monitoring approaches relevant to stock assessment of this species.


Introduction
Mesozooplankton (> 200 µm) such as copepods, have traditionally been considered the main prey of marine fish larvae (Munk and Kiørboe 1985;Pepin and Penney 1997). In the last decade, there has been growing recognition of the importance of smaller (< 200 µm) prey items such as phytoplankton, protozooplankton (PZP) and/or metazoan microzooplankton (MZP) (Montagnes et al. 2010;Denis et al. 2016;Bils et al. 2017) but their contribution to larval diets is still largely unknown. Laboratory studies have revealed that larval fish actively forage on PZP, such as Atlantic cod (Gadus morhua) (Hunt von Herbing and Gallager 2000), surgeonfish (Paracanthurus hepatus) (Nagano et al. 2000) and northern anchovy (Engraulis mordax) (Ohman et al. 1991). More recent laboratory studies on larvae of Pacific herring (Clupea pallasii) (Friedenberg et al. 2012) and Atlantic herring (Clupea harengus)  have reported that the presence of PZP in the diet can improve larval nutrition and survival. However, little is known at present on the relevance of PZP to larval fish diets in the field (Fukami et al. 1999;Figueiredo et al. 2005;Pepin and Dower 2007). The few studies conducted suggest that small prey can form a large portion of the larval diet (Figueiredo et al. 2007;Denis et al. 2016;Bils et al. 2017), especially during first feeding. The proportion of small prey decreases during ontogeny related to increasing mouth gape (Llopiz 2013 and references therein), larval length (Hufnagl and Peck 2011) and foraging ability (Munk 1992).
Assessing the contribution of PZP and MZP to larval fish diets is challenging as these small and often soft prey items are rapidly digested. Some species with long guts (e.g. clupeids) are also known to regurgitate prey during capture (Hay 1981). Moreover, gut content analyses of larval fish are typically performed on formalin-fixed material, a preservative unsuitable for PZP fixation. Therefore, indirect techniques, such as stable isotope analysis (SIA) offer a better chance to gain reliable estimates of the trophodynamic relationship of young larvae and their planktonic prey. The relative amount of stable carbon (δ 13 C) and nitrogen (δ 15 N) isotopes is commonly used to study food web dynamics, to detect the trophic position of an organism and to gain information on long term, integrated composition of diets (Peterson and Fry 1987;Vander Zanden and Rasmussen 1999;Post 2002;Boecklen et al. 2011). As prey items can be characterized by their specific stable isotope composition, habitat-specific differences in the feeding ecology of fish larvae (Laiz-Carrión et al. 2015;Costalago et al. 2016), changes in the dietary composition during ontogeny (Wells and Rooker 2009) as well as changing environmental conditions (Malzahn and Boersma 2009) can be revealed by SIA. These diet changes can be linked to larval nutritional condition and growth rates derived biochemically using the ratio of RNA to DNA (RNA/DNA) (Buckley et al. 2008).
For more than a century, it has been accepted that the survival and growth of early life stages form a bottleneck for marine fish recruitment (Cushing, 1975(Cushing, , 1990Iles and Sinclair 1982;Houde 2008). However, gaining a mechanistic understanding of the factors and processes driving fluctuations in year class success is still a major challenge in fisheries science. This is true even for species such as Atlantic herring, which were the focus of seminal works on recruitment variability and it is arguably the best studied small pelagic fish worldwide (Illes and Sinclair 1982;Peck et al. 2021). As other small pelagic fish, they have a boom-and-burst strategy, thus reacting quickly to environmental changes, and play a key role in the food web, linking plankton to higher trophic levels (e.g. adult fish or marine mammals). Particularly in the North Sea, herring populations have dramatically varied over the past decades, collapsing in 1980s and recovering in the 1990s. Current trends suggest that spawning stock biomass is relatively high, but recruitment continues to be relatively low (ICES 2021). In Autumn and Winter-Spawning herring in the North Sea (NSAS), recruitment appears to be driven by processes impacting the survival of overwintering larvae (Nash and Dickey-Collas 2005;Payne et al. 2009). Therefore, there is a special interest to disentangle the impact of different drivers to larval growth and survival in North Sea herring. Previous studies suggest that winter temperature, prey type and prey abundance play a vital role in year-class success (Fässler et al. 2011;Hufnagl and Peck 2011;Alvarez-Fernandez et al. 2015). The strength of any particular recruitment driver in North Sea herring may vary among spawning grounds since each component spawns at a unique time and location e.g. early/mid-September (Orkney-Shetland), late September (Buchan-Banks) and December/January (English Channel, Downs)   (Fig. 1). Larvae hatched at different seasons will target different prey items in relation to seasonal changes in the plankton community, and the contribution of PZP in larval diets would increase in winter when MZP abundances are lower.
The present study examined the contribution of PZP and MZP to larval herring diets using SIA. Furthermore, the relationship between prey field characteristics and biochemically based growth and condition of larvae was analyzed in two different years (2013 and 2014) at two spawning grounds and seasons in the North Sea (Buchan/Banks in autumn, and Downs in winter). Two hypotheses were tested: (1) larvae hatched later in the year (winter) rely more on PZP than on MZP to cover their energetic demands, and (2) larvae of different body sizes or that originate from different spawning grounds target different prey items. A better understanding of in situ prey availability, preferred prey items and potential implications for larval condition is needed to improve our knowledge on how environmental conditions impact early larval survival and the subsequent recruitment of herring in the North Sea.

Plankton and herring larvae sampling
Sampling of protozoo-, microzoo-and ichthyoplankton was conducted during the International Herring Larvae Survey (IHLS) on board RV Tridens. Samples were collected in September 2013 and 2014 on the Buchan and Banks spawning grounds (BB13 and BB14, Fig. 1) and in January 2014 on the Downs spawning ground in the English Channel (DO14). For certain analyses, Buchan and Banks spawning components are treated separately (Buc13, Buc14, Banks13 and Banks14; Fig. 1) to investigate potential differences between the autumn spawning components. Seawater samples were collected with Niskin bottles mounted on a CTD rosette sampler (Seabird SBE 911).
Samples for Chlorophyll a (Chl a) were taken at 15 m depth and samples for stable isotope analysis (SIA) of seston near the surface (max. 4 m depth) at 6 (BB13), 10 (BB14) and 8 (DO14) stations (Fig. 1). Up to 1000 mL of water was filtered through pre-combusted filters (Whatman GFC, 1.2-µm pore size). Duplicate filter samples were collected at Note CTD data of all stations were used for this plot, but plankton (seston, protozooplankton and microzooplankton in circles) and larval herring (crosses) were only sampled in selected stations for the present study. Note the temperature colour scaling (in °C) differs for Buchan/Banks and Downs 90 Page 4 of 16 each station for Chl a and seston and frozen at −80 °C until further analysis.
For PZP taxonomic analysis, a subsample from the Niskin bottles sampled near the surface was immediately stored in a 500-mL brown glass bottle and preserved with neutral Lugol's solution (2% final concentration).
A Gulf VII high-speed sampler (280 µm mesh size, 0.2 m nose cone opening) (Nash et al. 1998) and a PUP-net (52 µm mesh size) mounted thereon were used to capture herring larvae and microzooplankton (52-200 µm, MZP), respectively. The sampler was towed in double oblique hauls, at a speed of five knots, from the surface to 5 m above the seafloor. A sample of herring larvae was taken every 1/9 ICES rectangle (https:// www. ices. dk/ data/ maps/ Pages/ ICES-stati stical-recta ngles. aspx). The Gulf VII net was equipped with a flowmeter (Valeport) in the nosecone and a CTD (Seabird 911plus). Temperature and salinity were measured with the CTD mounted on the Gulf VII sampler at every station (n = 165, 145 and 96 stations in BB13, BB14 and DO14, respectively).
At the stations with PZP sampling, the MZP sample was split (Motodo Box plankton splitter) and half of the sample was separated into two size classes (52-100 µm and 100-200 µm) and filtered onto pre-combusted GFC filters for SIA. The other half of the sample was transferred to buffered 4% formalin for subsequent taxonomic analysis.
A minimum of 20 herring larvae were frozen at −80 °C onboard at 25 (BB13), 31 (BB14) and 20 (DO14) random stations where the total amount of larvae exceeded 100 individuals (based on onboard estimations) (Fig. 1). Frozen larvae were used for biochemical analysis (SIA and RNA/ DNA). Additional herring larvae caught with the Gulf VII sampler were preserved in 4% buffered formalin for length measurements. Larval length was measured to the nearest mm after preservation/freezing. No correction for shrinkage due to preservation and handling time was applied. A minimum of 20 larvae per station was sampled. The day/ night catch ratio for larvae < 10 mm is close to 1 (McGurk 1992). It increases with larval size but never exceeds 2.5 in larvae < 25 mm. As most larvae caught during the IHLS were < 12 mm (> 80% of total catches), day/night differences were not taken into account.

Plankton abundance
Chlorophyll was extracted with acetone (90%), incubated overnight and Chl a concentration was measured and calculated after Jeffrey and Humphrey (1975).
PZP seawater samples were analysed according to Utermöhl (1958). Due to the low PZP cell numbers, the samples were settled in a 100-mL sedimentation chamber (Hydro-Bios) for 48 h (HELCOM 2014). PZP cells were counted under an inverted microscope (Leica DMI 3000, 200 × with Moticam camera attached). The whole chamber plate was counted to avoid under-representation of less abundant PZP groups. PZP, comprising heterotrophic and mixotrophic ciliates and dinoflagellates, were identified to the lowest taxonomic level possible (Dodge and Hart-Jones 1982;Montagnes 1996;Olenina et al. 2006;Strüder-Kypke et al. 2006;Hoppenrath et al. 2009;Löder et al. 2012). Other protists and small-sized metazoans that occurred in water samples at very low abundance or at larger size (e.g. Foraminifera and copepod nauplii) were not counted and thus not included in the analysis. The PUP net samples (52-200 µm) were processed for PZP and MZP abundance with a FlowCam VS2 (Fluid Imaging technologies, Yarmouth, USA), an imaging particle analyzer, using a 4 × magnification. A subsample (dilution factor 10) of each sample was analyzed for BB13 and BB14, making sure that a minimum of 600 living organisms were counted (typically > 1000 organisms). The abundances (ind L −1 ) of the six most abundant planktonic groups were estimated: Tripos spp. (Dinophyceae), Dinophysis spp. (Dinophyceae), tintinnids (Ciliophora), copepods (including nauplii and copepodite stages) and mollusc larvae (Bivalvia and Gastropoda). In DO14, the whole PUP net sample was processed due to the relatively lower abundance of plankton in the samples. In DO14, only two groups (nauplii and small copepods/copepodites) were quantified, as other zooplankton organisms were very rare (Raudenkolb 2016). No copepod staging was applied.

Biochemical analysis
The GFC filters with seston and MZP were dried at 60 °C for 24 h prior to SIA and duplicate GFC filter samples were used for SIA of seston samples. SIA of MZP (52-100 µm and 100-200 µm) was performed on 1/4 of each filter to obtain the correct range of C and N. To estimate measurement error, two replicate quarters were analyzed separately, and mean values calculated.
For biochemical analysis, larvae between the lengths of 8 and 14 mm were chosen. This size range covers the life phase immediately after yolk absorption through the end of successful first feeding (Illing et al. , 2016. Due to the larger length-at-hatch, the size class S of larvae (8-9 mm) from DO14 might have contained yolk-sac larvae as well if they lost their yolk-sacs in the net. The 13 C and 15 N was analysed for 73 larvae (Table 1). Due to changes in C and N content with larval size (Kiørboe et al. 1987), SIA was analysed in small (S = 8-9 mm), medium (M = 10-12 mm) and large (L = 13-14 mm) length classes representing early, first, and advanced feeding larvae, respectively. At each station, larvae were pooled within these length classes to achieve the amount of C and N required for SIA. Larvae were combined from two (BB13 and BB14, respectively) or four (DO14) neighboring stations to obtain the necessary biomass for SIA measurements. Whole larvae were analysed without removing the guts since gut contents account for a maximum of 2% of total body weight (Pepin and Penney 2000). Larvae and filters were packed in tin capsules prior to SIA. All measurements of SIA ( 13 C and 15 N) were conducted at the Stable Isotope Facility at UC Davis, California, USA. The GFC filters (seston and MZP) were analyzed using an Elementar Vario EL Cube or Micro Cube elemental analyzer (Elementar Analysensysteme GmbH, Hanau, Germany) interfaced to a PDZ Europa 20-20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK). The seston filters contained a mixture of phytoplankton, detritus and PZP. Samples containing < 100 µg C or < 20 µg N were excluded from further analysis according to the sensitivity ranges given by the Stable Isotope Facility at UC Davis, California, USA.
The 13 C and 15 N isotopes of the herring larvae were analysed using a PDZ Europa ANCA-GSL elemental analyser interfaced to a PDZ Europa 20-20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK). δ 13 C values of the larvae were corrected for lipid content as suggested by Post et al. (2007) for a C-N ratio > 4: The stable isotope ratios were expressed as conventional δ notation where δX = 15 N or 13 C, R = 15 N: 14 N or 13 C: 12 C and the relative international standards of Vienna PeeDee Belemnite for nitrogen (± 0.3‰ SD) and Air for carbon (± 0.2‰ SD).
Some of the remaining frozen larvae were used to estimate nutritional condition (RNA/DNA). In total, 117 (BB13), 95 (BB14) and 171 (DO14) herring larvae were classified as intact and larval length was measured to the nearest mm with a stereomicroscope (Leica MZ16). Afterwards, larvae were freeze dried at −50 °C (Christ Alpha 1, 4 LSC) and dry mass (DM, µg) was determined using a microbalance (Sartorius Genus SE2, ± 0.1 µg). The RNA/ DNA was determined on single larvae following a modified protocol of Caldarone et al. (2001) that uses ethidium bromide as binding fluorescent dye and restriction enzymes for RNA and DNA. The standardized RNA/DNA (sRD, Caldarone et al. 2006) was estimated using a factor of 2.4 based on the slopes of the RNA and DNA calibration series.
The instantaneous growth rate ( Gi , d −1 ) was calculated from sRD and in situ water temperature (T, °C) using the following equation (Buckley et al., 2008):

Data analyses
All statistical and graphical analyses were performed using the R software (V 3.1.0 R core team 2014). Sea surface temperature (SST, 5 m depth) was mapped using inverse distance interpolation in the RGeostats package (Renard et al. 2016). Thermal stratification was estimated as the temperature difference of min. 1.5 °C between SST and bottom temperature.
Relationships between environmental variables (SST, salinity) and the abundance of herring larvae and PZP were explored using Pearson correlation. Similarly, this test was also used to explore relationships between larval herring abundance and prey fields (abundances of ciliates, dinoflagellates, copepod nauplii and mollusc larvae).
The δ 15 N of the seston samples (δ 15 N Seston ) was used as a baseline to estimate the trophic position (TP) in a one-source food web model (Hobson and Welch 1992). A trophic fractionation factor (TFF) of 2.2‰ (invertebrates) and 3.4‰ (vertebrates) between trophic levels for δ 15 N was assumed  (Minagawa and Wada 1984;Peterson and Fry 1987;Post 2002): To test for significant differences in growth rate and sRD of the larvae and δ 15 N and δ 13 C of larvae and plankton components, one-way ANOVA's were performed within as well as among larval length classes, spawning grounds and spawning seasons followed by a Tukeys post hoc test. Some of the size classes did not meet the requirements for ANOVA, so a non-parametric test (Kruskal Wallis, followed by a Dunn's Test) was used in these cases.
The sRD was normalized by performing a square-root transformation. For each spawning ground, a one-way ANOVA (followed by a Tukey post hoc test) was performed on the transformed sRD data and a multivariate linear regression model (GLM) was used to reveal potential effects of larval length and/or environmental variables (SST and salinity). Larvae were grouped into the three length classes used for SIA (S, M and L) (see previous section). Samples collected at BB13 and BB14 were from two different spawning grounds, Buchan in the north (Buc13 and Buc14) and Banks in the south (Banks13 and Banks14) ( Figure S1). Therefore, the areas are treated as separated variables for statistical analyses of nutritional condition, growth rate and stable isotope composition.

Environmental conditions and larval prey abundances
In BB13, sea surface temperature (SST) ranged from 11.6 °C to 14.7 °C ( Fig. 1C) with colder waters in the north compared to the south. The same north-south pattern was observed in BB14 (Fig. 1D), but SST was higher and ranged from 13.0 to 15.2 °C. In DO14, winter SST was between 6.9 and 10.9 °C (Fig. 1B). A salinity gradient occurred in the English Channel with higher salinities in the western and lower salinities in the eastern part of the sampling grid (29.2-35.2). In BB13 and BB14, the surface salinity ranged from 34.0 to 35.1 and 34.0 to 34.9, respectively, with the highest salinities in the northern part of the sampling area. In BB13 and BB14, the water column was thermally stratified (at 58 and 81 stations, respectively) in contrast to the winter conditions during DO14, where the water column was well mixed in all stations (not shown). The Chl a concentration was more variable in BB13 and ranged from 1.39 to 8.95 µg L −1 compared to BB14 (range 2.2 to 4.86 µg L −1 ) ( Table 2).
In contrast to autumn, the PZP abundance during wintertime (DO14) was one order of magnitude lower (with a maximum abundance of 1750 ind L −1 ), but the community still comprised mainly of small cells (< 50 µm). Some taxa common in BB13 and BB14 were rare in DO14, such as P. pelagica, all species of Tripos and Dinophysis sp.. The aloricate ciliates of the genus Strombidium (range 10 ind L −1 to 220 ind L −1 ) and Strobilidium (range 10 ind L −1 to 170 ind L −1 ) dominated the ciliate community. During winter, tintinnids were the third most abundant ciliate taxon contributing up to 25% to total ciliate abundance compared to 10% in autumn. Total PZP abundance (as well as ciliate and dinoflagellate abundance) was positively correlated with SST and chlorophyll a (Pearson correlation, p < 0.02) when the data from different spawning grounds were pooled and analyzed. However, when the different areas were analyzed separately, there were not significant relationships (Pearson correlation, p > 0.1).The abundance of copepod nauplii in BB13 was higher in Buc13 than Banks13 (max.30 vs. max. 8 ind L −1 , respectively), whereas the opposite pattern was found in BB14 (max. 13 vs. max. 23 ind L −1 in Buc14 and Banks14, respectively). Copepod (copepodites and adults) abundance ranged from nearly 0 to 1 ind L −1 in BB13 and   BB14. Mollusc larvae were rare in BB13 (< 1 ind L −1 ) but reached up to 536 ind L −1 in BB14. During wintertime (DO14), copepod nauplii were present at all stations (< 8 ind L −1 ), and copepodites were also present but at lower abundance (< 1 ind L −1 ) ( Table 2).

Herring larvae abundance and distribution
In BB13 and BB14, larval herring were mostly found in the northern (Buchan) and southern (Banks) regions of the survey grid ( Figure S1). Peak abundance was higher in BB13 (270 ind m −3 ) than BB14 (13 ind m −3 ). In both years, most larvae were 7 to 9 mm in length ( Figure S2). In DO14, the maximum abundance of larvae was 50 ind m −3 and consisted mostly of 11 to 12 mm individuals. It must be noted that the length of the frozen larvae used for the analysis presented here cannot be directly compared to the counts conducted for the ICES length-frequency calculations, where the length of formalin-preserved larvae was used.
Larval abundance was unrelated to the station-specific abundance of ciliates, dinoflagellates, copepod nauplii, and mollusc larvae at each spawning ground and overall (Pearson correlation, p > 0.23). Using the complete dataset, no consistent relationship between larval abundance and environmental variables was observed either. SST was only weakly related to larval abundance in BB14 (Pearson correlation, σ = −0.19, p = 0.02), but not in the other spawning grounds (Pearson correlation, p > 0.1). Salinity was unrelated to larval abundance in all cruises (Pearson correlation, p > 0.5).

Herring larvae trophic position
Seston samples had an average δ 15 N of < 7‰ on every survey (Fig. 2). A significant difference in δ 13 C between the plankton size classes (50-100 µm and 100-200 µm) was found only for Buc14. Thus, the size classes were pooled for the analysis and are further referred to as MZP (Table 1). MZP was significantly more enriched in δ 13 C (p < 0.001) in Buc13 and DO14 compared to Banks13 (Fig. 2).
The δ 15 N of the larvae showed no statistical differences between the Buc13 and Banks13 for small and large larvae (ANOVA, p > 0.05). In length class M, lower δ 15 N was observed in Buc14 compared to Banks14 (ANOVA, p = 0.049). In DO14, δ 15 N of all length classes were higher and differed significantly from the values in BB13 and BB14 (ANOVA, Tukeys post-hoc p < 0.001, Length class S Kruskal-Wallis, Dunn's test p < 0.01). Note that δ 15 N estimates were not available for Banks13 size M, Buc14 size L and Banks14 size L (Table 2).
In Buc13, small larvae were more depleted in δ 13 C than large larvae (ANOVA, p < 0.001) and showed higher δ 13 C than larvae at Banks13 for all size classes (ANOVA p < 0.001) (note: no larvae were collected for length class M). For δ 13 C, significant differences in DO14 compared to BB13 and BB14 were found between larvae of length class M (ANOVA followed by Tukey post-hoc, p < 0.001, p = 0.02, respectively). Note that δ 13 C estimates were not available for Banks13 size M, Buc14 size L and Banks14 size L (Table 2). MZP was less than one trophic level above baseline (Fig. 2) and similar for all surveys with a trophic position (TP) of 1.27 (BB13), 1.19 (BB14) and 1.17 (DO14). Larvae of all length classes in BB13 showed a TP between 2.08 and 2.28, independent from the spawning ground.
The trophic enrichment (∆ 15 N) relative to the SI of MZP was in the range of the estimated enrichment of 3.4‰ between two trophic levels (3.07-3.39‰) except for the largest larvae in Banks13 (2.74‰). In BB14, ∆ 15 N relative to MZP was higher and reached values exceeding one trophic level (up to 4.39‰ in Banks14). The larvae in DO14 showed an enrichment of > 5‰ relative to MZP (Fig. 2).
The relation of the total amount of carbon to nitrogen in the larvae decreased significantly with larval dry weight in BB13 (C:N = 0.0003 dry weight + 3.9242) (Fig. 3). Larvae of BB14 and DO14 showed no significant relationship between C:N ratio and larval weight. Larvae of length class

Herring larvae nutritional condition and growth rate
Larval nutritional condition (sRD) did not significantly differ between the length classes within each spawning area (Fig. 4) (ANOVA, p > 0.05), except for larvae of length class M in Banks14, which showed a significant lower sRD than the other larvae in BB14 (Tukey post hoc, p = 0.03). The lowest sRD was found in DO14 with 0.42 in length class M, the highest in Buc14 with 6.39 in length class M (Fig. 4). Overall, temperature had a strong positive effect on larval sRD (GLM, p < 0.001) (not shown), but no correlation was found within each single survey. Salinity was not related to sRD (GLM, p > 0.05).
Regarding interannual and seasonal differences, growth (G i ) results are reported for a better comparison as they account for thermal differences. All larvae in the three surveys showed positive G i (> 0.02 d −1 ). For size classes S and L, larvae in DO14 showed a significant lower G i than in Buc13 (Tukey post hoc, p = 0.044) with maximum G i of 0.4 d −1 . For size class M, Banks14 and DO14 had a significantly growth rate compared to the other areas and years (Tables S2 and S3).

Discussion
In environments with strong seasonal differences in biological productivity, periods of low primary production and low metazoan plankton abundance, protozooplankton is considered to be an important prey item for marine fish larvae, especially during first feeding (Montagnes et al. 2010). However, this assumption has seldomly been tested. Here, we used North Sea herring larval surveys to test whether differences in the availability of PZP and MZP can alter larval fish diets and nutritional condition. PZP and MZP communities were similar in both seasons in terms of diversity and isotopic signature (δ 15 N), although the abundance was lower during wintertime. Our results do not suggest that PZP has a higher contribution than MZP to the diets of winter larvae, and, in fact, δ 15 N was higher in these compared to autumn larvae, probably related to the maternal signature. (10-12 mm) and L (13-14 mm). Boxplot colour in BB13 and BB14 refers to subareas: white for Buchan (Buc13 and Buc14) and grey for Banks (Banks13 and Banks14). The number of larvae used for the analysis is shown. Data with n < 3 were not used for the analysis Surprisingly, winter larvae (DO14) grew at the same pace as autumn larvae (BB14), based on the biochemically-derived growth rate (from RNA/DNA estimates). There were no significant changes in trophic position through ontogeny (larvae 8 to 14 mm), but there was a tendency that larvae from different size classes fed on different MZP components from the same trophic level (similar δ 15 N but different δ 13 C). Altogether, these results suggest that herring larvae feed on a diverse array of PZP and MZP prey items during both seasons. The limitations of the bulk SIA approach used in this study are discussed in the light of recent multi-proxy studies on trophic interactions.
In the following, our two hypotheses will be discussed: (i) larvae hatched later in the year (winter) rely more on PZP than MZP to cover their energetic demands, and (ii) larvae of different body sizes or hatched at different spawning grounds target different prey items.

Larvae spawned in winter rely more on PZP than MZP to cover their energetic demands
The recorded environmental conditions on the spawning grounds studied here were comparable to those reported in long-term time series data from North Sea monitoring sites and CPR data (Continuous Plankton Recorder) (www. wgze. net, Bresnan et al. 2015). Bottom temperature at DO14 and BB14 was about one degree higher than in previous years (van Damme and Bakker 2014) which is in line with positive temperature anomalies observed since mid-1990s (O'Brien et al. 2013).
In wintertime (DO14), productivity was lower compared to autumn (BB13, BB14) with Chl a concentrations < 1.5 µg Chl a L −1 (www. wgze. net, this study) and PZP abundances < 1750 ind L −1 . Both periods displayed a similar composition with athecate dinoflagellates dominating the community as opposed to ciliates observed in other studies in the region (Widdicombe et al. 2010;Löder et al. 2012;Bils et al. 2019). The total abundance of PZP and nauplii found here was in the range of those recorded in long-term sampling campaigns in the English Channel (Widdicombe et al. 2010) and at a monitoring site in the eastern North Sea (Wesche et al. 2007). This homogenous composition in autumn and winter was reflected in the isotopic signatures of the seston and MZP samples. δ 15 N was similar in seston and MZP during all three surveys. However, δ 13 C of the MZP was elevated in Buc13 and DO14 compared to the other samples. One possible explanation for this elevated δ 13 C was that the organisms of the MZP fraction in those areas possessed higher lipid storage (and thus higher nutritional value) (Waite et al. 2007).
Despite the absence of seasonal differences in the seston and MZP δ 15 N, winter larvae (DO14) showed elevated δ 15 N compared to autumn-spawned larvae. This pattern was opposite to what we initially expected: that δ 15 N in winter larvae would be lower than in autumn. The PZP community in this study was dominated by omnivorous, phagotrophic or bacterivorous organisms (e.g. Heterocapsa sp., Gymnodiniales, Strombidium spp.). Thus, no significant contribution of carnivorous unicellular organisms was observed, which could have lead to higher trophic position of the larvae compared to when primarily herbivorous copepods were the main food source (Pepin and Dower 2007). One explanation for this higher δ 15 N in the winter larvae may be related to a higher contribution of the maternal signature. Adult herring tend to have higher δ 15 N in coastal and southern areas of the North Sea (Jansen et al. 2012), where the winter spawning grounds are located. Winter larvae also spawn with a larger size (see "Herring prey items depend on larval size and region"). Thus, the contribution of maternal signals will be stronger in these small and medium sizes (9-12 mm). The results of the present study do not support our initial hypothesis that winter larvae rely more on PZP than MZP.
Our results suggest a disconnection between the δ 15 N and δ 13 C of seston and MZP and the potential prey items consumed by the larvae. One potential explanation for the δ 15 N enrichment in winter-spawned larvae that is not reflected in the seston and MZP δ 15 N might be related to the fact that isotope enrichment processes happen to different degrees for heterotrophic protists (Gutiérrez-Rodríguez et al. 2014). These authors showed no or only little trophic fractionation between heterotrophic protists and their autotrophic prey occurs, in contrast to metazoan organisms. Thus, the 2.2‰ enrichment usually observed for invertebrates and 3.4‰ enrichment for vertebrates may not necessarily be applied for consumers feeding on protozoa. Overall, the δ 13 C of MZP in our study are above average for pelagic plankton organisms (Kürten et al. 2013). This might be partially explained by higher lipid content of the MZP causing an elevation of δ 13 C (Waite et al. 2007) or due to a higher contribution of benthic components in turbulent environments (Marconi et al. 2011). Regarding the latter, wind-driven mixing may have suspended benthic detrital material and sediments into the plankton that were sampled in the MZP net, as some benthic cells were visually observed in the samples (Justus van Beusekom and Albert Calbet, pers. comm).
Based on our results, PZP seemed to contribute to the diets of herring larvae. Whether the contribution of PZP prey is of higher importance to the larvae in winter compared to autumn (when metazoan plankton is more abundant) remains unknown. New available methods may help to shed light on the link between PZP and fish larvae, as the methods used in this study showed some limitations (see "Methodological assessment"). Given the importance of prey (< 300 µm) in larval herring diets (Alvarez-Fernandez et al. 2015;Dudeck et al. 2021), further studies on this trophic link are needed, as well as considerations to include PZP and MZP monitoring into routine programs to improve current fisheries assessment processes. For example, if PZP and MZP assessments suggest that certain taxa or size classes are directly related to larval growth, they could serve as indicators of good habitat suitability and larval fish performance (e.g. Batten et al. 2016).

Herring prey items depend on larval size and region
The present study showed no differences in δ 15 N in herring at different larval sizes (8 to 14 mm) collected during autumn and winter in the North Sea. This suggests that larvae in all size classes preyed on the same trophic level, and that larvae consumed a wide range of prey items, not being specifically more selective at larger body sizes (Malzahn and Boersma 2009;Wells and Rooker 2009). The δ 15 N values measured for BB13 and BB14 are in accordance with experimentally derived data on post-feeding herring larvae (Aberle and Malzahn 2007) but lower than mean δ 15 N known for adult herring in the North Sea (e.g. Jansen et al. 2012). Therefore, the δ 15 N of herring larvae from the Buchan-Banks area (Autumn) clearly provided a first/ advanced feeding signal, as opposed to yolk-sac larvae that usually reflect the maternal signal.
Differences in δ 13 C between small larvae (8-9 mm) and larger larvae (13-14 mm) within the Buchan and Banks spawning grounds make it plausible that the larvae used different primary carbon sources (difference in δ 13 C) from the same trophic level (no differences in δ 15 N). Trophic enrichment (∆ 15 N) in BB13 did not show an increase with larval size, thus suggesting that the larvae of all length classes preferred prey from the MZP size fraction rather than from the seston. Larvae from Banks14 achieved higher ∆ 15 N enrichment that exceeded the estimated enrichment of 3.4‰ per trophic level thus pointing at a consumption of prey items with a higher trophic position than MZP. Thus, the bulk plankton signals may not specifically reflect the prey spectrum these larvae fed on (see Sect. "Methodological assessment"). Larvae spawned in winter (DO14) achieved unexpected high δ 15 N, as discussed in the previous section, and a higher enrichment than normally expected for one trophic level. Nevertheless, the size classes within DO14 did not significantly differ from each other in isotopic composition, indicating a similar prey field of larvae between 8 and 14 mm.
Our initial hypothesis was that there will be differences in prey consumed (derived from SIA analysis) throughout ontogeny in both spawning grounds. Such shifts in larval diet during ontogeny are known and have been previously revealed by SIA of various marine fish species in the Gulf of Mexico (Wells and Rooker 2009), such as during the settlement process in red drum (Sciaenops ocellatus) (Herzka et al. 2002). For NSAS herring in Downs, Denis et al. (2016) reported a shift in diet at a length of ~ 12 mm using scanning electron microscopy (detecting microplanktonic organisms in larval guts, but without accounting for naked PZP). However, our study did not confirm those SIA results. This might be related to limitations of the current bulk SIA analysis pooling several larvae (see Sect. "Methodological assessment"), as C:N was inversely related to dry weight for larvae in BB13. These results suggest that the larvae were feeding on qualitatively higher food sources as they grow, as shown in other studies in larval herring. For example, Ehrlich (1974) reported that C:N in exogenously feeding herring larvae decreased with development from 4.05 at the end of the yolk-sac stage to 3.63 for 21-mm larvae.
One important aspect that must be considered while comparing different spawning grounds is that larval length at hatch varies. Both autumn-and winter-spawning herring start oocyte development at the same time (March/April). Oocytes will continue growing until they are layed, leading to larger eggs (and thus larvae) in winter compared to autumn (Damme et al. 2009). Downs larvae usually hatch at a length of 9 mm (Heath et al. 1997) and Buchan and Banks at 7 and 8 mm, respectively (Blaxter and Hempel 1963). In addition, growth and development of larval fish is slower at colder temperatures (Pepin 1991). Larval sizes were larger in September 2014 than 2013 during the same sampling period, indicating that spawning had occurred earlier that year, potentially induced by higher overall temperatures (Blaxter and Hempel 1963). However, higher temperatures and higher MZP abundance in autumn 2014 compared to 2013 did not promote a positive effect on growth and condition nor did it alter δ 15 N or δ 13 C, as larval growth (mean G i ) for larvae between 10 and 14 mm was significantly higher in autumn 2013 compared to 2014. The abundance of small copepods and copepodite stages was lower in BB14 compared to BB13 and DO14, below the suggested minimum prey threshold (Munk and Kiørboe 1985;Figueiredo et al. 2005;Hufnagl and Peck 2011) to sustain larval growth. Growth rates during winter were > 0.02 d −1 , well above the reported growth rates in January 2015 (< 0 d −1 , Denis et al. 2017).
Even though SIA (δ 15 N and δ 13 C) of herring larvae and bulk plankton samples revealed no differences in trophic level between larvae from 8 to 14 mm within each spawning ground, differences in larval growth and condition were observed in autumn-spawned larvae suggesting that feeding on a higher trophic level is not directly linked to higher growth rates or better nutritional condition. The ∆ 15 N led to the assumption that the larvae hatched in autumn are covering their nutritional needs mostly with MZP prey. Winterspawned larvae in the English Channel exhibited δ 15 N > 1‰ higher than autumn-spawned larvae. However, the higher trophic position of larvae could not directly be related to the SIA of prey items.

Methodological assessment
The present study gives some indications about different feeding habits among larval sizes and spawning grounds, but also shows limitations of the bulk SIA approach. Bulk samples of plankton size classes of potential prey do not account for the high variability of trophic modes (relative contribution of autotrophic and heterotrophic fractions and non-living detrital material) and species-specific stable isotope signals of prey (Kürten et al. 2013;Schoo et al. 2018). Recent SIA approaches proposed a thorough sorting of specific plankton groups/species to identify SIA of prey instead of using bulk seston or rough size-fractionated plankton fractions (Costalago et al. 2020) and the use of multi-proxy approaches to gain a broader understanding of trophic dynamics is stressed (Amundsen and Sanchez-Hernandez 2019;Bachiller et al. 2020). Besides SIA, fatty acid analyses could complement food web analyses using a dual-tracer approach (Hielscher et al. 2015;Schoo et al. 2018), while DNA barcoding of gut contents is suggested as an alternative to overcome the problem to identify prey in the digestive tract (Roslin et al. 2016). This approach was applied successfully, for example, in larvae of European eel (Anguilla anguilla) (Riemann et al. 2010;Ayala et al. 2018). This and other genomic techniques applied in recent years to discover PZP prey in gut contents also bear caveats, ranging from varying digestion times of different PZP species to bias by protist symbionts of the predator (Stoecker and Pierson 2019 and references therein). As reported in Bachiller et al. (2020), combined approaches using SIA and molecular gut content analysis have a high potential for detailed food web analysis with a high species-specific resolution and are thus the next step forward to characterize dietary composition, feeding preferences and contribution of small prey to the diets of fish larvae.
Besides the bulk SIA, it is worth noting other methodological weaknesses that could have impacted the results. The number of SIA samples per larval size class was unbalanced (between 3 and 12 per size class, Table 1), due to difference in the number of larval collected within each size class during the surveys. This unbalanced sampling and the pooling of several larvae per sample (up to 7) could have impacted our results in terms of the changes in diet through ontogeny and comparisons across spawning grounds. It is important to note that the proxies used here generally reflect larval foraging and/or growth over a few days prior to sampling (e.g. 2-3 days for RNA/DNA, Peck et al. 2015), so that it is not unexpected to see no relationship with the in situ plankton abundances. Additionally, seasonal comparisons need to consider that the geographic location of the spawning grounds change, not allowing for a direct comparison in the exact same location. Finally, Downs was only sampled once, which prevented from exploring interannual variability for this spawning ground.

Conclusion and outlook
Despite the limitations of the bulk SIA approach, our results suggest that larval herring in the North Sea rely on a diverse diet based on PZP and MZP in winter and autumn. Given the current concerns about relatively low recruitment success in North Sea herring and the increasing contribution of the Downs component, further studies using multi-proxy trophic approaches that allow a more detailed analysis of prey selectivity at a higher speciesspecific resolution are needed. Such information combined with the implementation of PZP and MZP during routine larval surveys for fisheries assessments (e.g. IHLS), can provide important information to understand interannual and seasonal changes in prey fields and how these relate to larval growth and survival. For example, diatoms and loricate ciliates explained 75% and 45% of the interannual variability in recruitment of Pacific herring in Prince William Sound (Batten et al. 2016). Combining this knowledge from statistical studies with in situ field work to explore trophic interactions can help develop indicators of prey abundance from relevant taxa. Such indicators and knowledge can later support the stock assessment process in the calculations of the biological reference points used in the assessment models. Moreover, the fast development of climate forecasting brings hope that the application of these indicators could be soon implemented at time-scales relevant to provide advice to stock assessment and ecosystem-based management (Hobday et al. 2018).