Metabarcoding analysis suggests that ﬂexible food web interactions in the eukaryotic plankton community are more common than speciﬁc predator–prey relationships at Helgoland Roads, North Sea

Metabarcoding analysis suggests that ﬂexible food web interactions in the eukaryotic plankton community are more common than speciﬁc predator–prey relationships at Helgoland Roads, Sea. Various ﬁeld studies on plankton dynamics have broadened our understanding of seasonal succession patterns. Additionally, laboratory experiments have described consumers ranging from generalists to selective grazers. While both approaches can give us a good understanding of the ecosystem and its dynamics, drawbacks in identiﬁcation and a limited coverage of the ecosystem have left open questions on the generality of previous results. Using an integrative approach, we investigated water samples taken at Helgoland Roads by metabarcoding to describe seasonal succession patterns of the whole plankton community. By use of network analysis, we also tried to identify predator–prey dynamics. Our data set depicted the strong seasonality typically found for temperate waters. Despite a stable background community surviv-ing strong ﬂuctuations, small and abrupt changes, such as pronounced blooms and random appearance of autotrophs, cause seasons to be quite different in an inter-year comparison. Main consumers were copepods, ciliates, and dinoﬂagellates, of which the latter were most abundant. Furthermore, our results suggest that zooplankton predators might favour speciﬁc prey during certain time periods but seem to be quite opportunistic otherwise throughout the year. The variability and potential for many different relationships in the plankton community might be an indicator of resilience in the system.


Introduction
The immense diversity and size spectrum of marine eukaryotic plankton makes it difficult to capture and study the whole community at once. As a result, our understanding of marine ecosystems is somewhat fragmentary, as field studies out of methodological necessity focus typically on individual community compartments, either in temporally explicit one-point timerepeated measurements or spatially explicit (one-time many points) settings. The findings from these field studies, or the questions arising from them, are typically then addressed using laboratory, mesocosm, or whole-field experiments. The use of frequent monitoring at Long-Term Ecological Research (LTER) sites is another approach to broaden the view on seasonal and longerterm scales.
The unique observation programme Helgoland Roads studies the long-term development of abiotic factors such as temperature, and the resulting dynamics of the planktonic community, at the Helgoland Roads station in the North Sea . This programme has been boosted by additional studies focusing on different, specific, plankton groups to provide a more detailed view on their dynamics (e.g. Wiltshire and Dürselen, 2004;Medlin et al., 2006;Sapp et al., 2007;Knefelkamp, 2009;Metfies et al., 2010;Löder et al., 2011Löder et al., , 2012Schlüter et al., 2012;Yang et al., 2014;Boersma et al., 2015;Lucas et al., 2015;Wiltshire et al., 2015;Scharfe and Wiltshire, 2019). These studies helped to improve the understanding of the seasonal dynamics of the different food web compartments at Helgoland. The abundance of picoeukaryotes, for instance, was found to be the highest in spring and summer and rapid changes could be observed throughout the year at Helgoland (Knefelkamp, 2009). It was also shown that the contribution of cryptophytes to the picoplankton was the highest in winter and spring (Metfies et al., 2010). Autotrophic plankton such as diatoms is known to be highly abundant during spring, when environmental conditions such as temperature and light availability increase (Mieruch et al., 2010). While heterotrophic dinoflagellates were found to be the most important contributors to biomass in general, mixotrophic dinoflagellates and ciliates can also significantly contribute to the total planktonic biomass at certain times of the year (Löder et al., 2012). Generally, microzooplankton might exert a stronger topdown control on phytoplankton at Helgoland than mesozooplankton (Löder et al., 2011).
Whereas the above-described specific studies on plankton trophodynamics have certainly increased our knowledge of relatively short-term patterns, or of single components of the food web, an overall assessment of the complete planktonic component is still lacking. Furthermore, by focussing on conspicuous or short-term food web interactions these might be over-interpreted, if they are assumed to be a regular phenomenon. New technology and analytical approaches might fill this knowledge gap by providing information on the complete plankton community as well as on interactions between different food web components, such as between nano-and picoplankton, and micro-and macroplankton. For instance, metabarcoding approaches have been used to study planktonic microbial communities around Helgoland, allowing the identification of new succession patterns in bacterioplankton throughout the year (Chafee et al., 2018), and in small eukaryotic plankton during spring .
Ideally, metabarcoding studies should be used to study the whole ecosystem at once, integrating all components. These integrated approaches are, however, rare (but see Abad et al., 2016Abad et al., , 2017. Instead, most metabarcoding studies either focus on smaller components of the community, by studying water samples, after removing larger mesozooplankton (Hernández-Ruiz et al., 2018;Rachik et al., 2018;Giner et al., 2019;Sprong et al., 2020) or they study the larger zooplankton, typically using net samples (Lindeque et al., 2013;Sommer et al., 2017;Bucklin et al., 2019;Blanco-Bercial, 2020). Here, we aimed to integrate these approaches by investigating the whole planktonic community using metabarcoding at once.
With this method, knowledge of zooplankton biodiversity, which is probably much higher than known today, and the functional role of different zooplankton species in the planktonic food web can be further enhanced. The links in the planktonic food web, identified by metabarcoding, can be visualized by conducting network analysis (Kurtz et al., 2015). Hereby, interactions or associations (for example, in terms of co-occurrence) are shown by edges (also called links) that connect different nodes (also called vertices). By using these kinds of techniques, we can break down complex community structures into distinct clusters at different times. These clusters can then be used to obtain new insights into the relationships within the planktonic community throughout the food web, and to discover potentially new or to verify previously observed predator-prey relationships.
Hence, for the first time, we conducted a metabarcoding study over a 3-year period to provide a comprehensive assessment of the annual succession of species constituting the planktonic food web presented in 1 L of water. In addition to identifying plankton community diversity as a whole, we suggest that metabarcoding analysis of natural seawater could provide unique insights into potential predator-prey relationships within the planktonic food web. Our aims were (i) to identify plankton communities that occurred from 2016 to 2019 and their succession using metabarcoding analysis and (ii) to identify predator-prey dynamics with regard to zoo-and phytoplankton. We used information on previously observed predator-prey pairings to check for co-occurrences of consumers/predators (in the following only referred to as predators) and prey in the metabarcoding data set. Furthermore, conducting of network analyses gave us a unique possibility to investigate associations and to corroborate potential relationships in the plankton community.

Study site and sampling
Work daily water surface samples (depth: 1 m) were taken with a bucket at the Helgoland Roads LTER sampling site (54 11 0 N, 7 54 0 E) from mid-March 2016 to mid-March 2019. The LTER site is situated between the main island and the dune island of Helgoland ( Figure 1). Due to the strong tidal currents, the surface samples are representative of the complete water column at the station (Hickel, 1998;Wiltshire et al., 2015). Depending on the tides, the well-mixed water column can fluctuate between 10 m depth (Callies and Scharfe, 2015). Measuring of Secchi depth and temperature was conducted directly at the sampling site. Other parameters were measured in the laboratory according to the LTER protocols (Hickel et al., 1993;Wiltshire et al., 2008Wiltshire et al., , 2010. These include salinity and nutrients such as silicate, phosphate, and inorganic nitrogen (nitrate, nitrite, and ammonium; Grasshoff, 1976) and chlorophyll a. Daily observations of sunshine duration in hours were downloaded from the Deutscher Wetterdienst, Climate Data Centre (DWD Climate Data Center (CDC), 2019); and the seasons were defined using the meteorological calendar: Spring ¼ March to May, Summer ¼ June to August, Autumn ¼ September to November, Winter ¼ December to February.

Molecular analysis
We combined samples from three different sampling periods for metabarcoding analysis. The sampling protocols on sampling frequency, filtration, and DNA extraction steps of the different sampling periods were not identical. However, each sample was taken with a bucket and 1 L of seawater was filtered. The first sampling period from 15 March 2016 to 31 May 2016 included work daily sampling and a sequentially filtration using 10 mm polycarbonate filters, 3 mm PC filters, and 0.2 mM polyvinylidene fluoride filters (Millipore, Schwalbach, Germany; Teeling et al., 2016;Käse et al., 2020). The samples from the other two sampling periods (from June 2016 to March 2019) were filtered onto 0.45 mm nylon filters (Whatman, 47 mm). Samples from June to October 2016 were taken infrequently (in total six samples; Sprong et al., 2020). From December 2016 until 14 March 2019, the samples were analysed twice a week, with additional samples from mid-May to the end of July 2018 (three samples per week). All samples were stored at À20 C until DNA isolation.
General methods on DNA extraction, MiSeq TM Illumina sequencing, and data processing have been described elsewhere (Käse et al., 2021).
In short, two different protocols were used for DNA extraction. DNA extraction for the 10 mm, 3 mm filters of the sampling period in spring 2016 and all 0.45 mm filters from June 2016 to March 2019 was carried out with the Macherey-Nagel NucleoSpin V R Plant II Kit. DNA extraction from 0.2 mm filters from spring 2016 was conducted as described previously by Sapp et al. (2007). It should be mentioned that DNA of multi-celled zooplankton or other big organisms is o based on DNA that occurred in the 1 l water sample. This can include DNA sticking to particle surfaces or faecal pellets as well as free DNA. No additional samples of these big organisms were included and macroalgae, copepods, or other mesozooplankton, that were visible to the naked eye, were not present during the DNA extraction process.
Sequentially filtered samples were pooled accordingly to achieve one sample per sampling date. A fragment (V4 region) of the 18S ribosomal (r) DNA was amplified using KAPA HiFi HotStartReadyMix (Kapa Biosystems, Inc., Roche, Germany) and the following primer set developed by Fadeev et al. (2018): 528iF (GCG GTA ATT CCA GCT CCA A) and 964iR (ACTTT CGT TCT TGA TYR R). About 43 million 2 Â 300 bp paired-end sequences were produced using an Illumina MiSeq TM sequencer (Illumina, San Diego, CA, USA). After data processing, 21 million sequences remained and were clustered using swarm (version 2.2.2, Mahé et al., 2014Mahé et al., , 2015. The Protist Ribosomal Reference database (PR2), version 4.11.1 (Guillou et al., 2013) was used as reference and all names are based on the taxonomy as it is given by the database. Sequence data for this study have been deposited in the European Nucleotide Archive (ENA) at European Bioinformatics Institute (EMBL-EBI) under accession number PRJEB37135, using the data brokerage service of the German Federation for Biological Data (GFBio, Diepenbroek et al., 2014), in compliance with the Minimal Information about any (X) Sequence (MIxS) standard (Yilmaz et al., 2011). Details of our pipeline are available on GitHub (https://github.com/PyoneerO/ qzip) and the full table of operational taxonomic units (OTU280 samples with 59.284 OTUs) was archived in PANGAEA (https:// doi.pangaea.de/10.1594/PANGAEA.921026).
A threshold of 0.001% (of total reads) was applied to the full data set, resulting in a data set of 2790 OTUs, which was used for all further analyses. Identification up to genus level was accepted. Where necessary, identification on species level was verified with BLAST. For taxa that could not be further identified, the previous taxonomic level was assigned; these objects were indicated with an additional term (e.g. unclassified) and they were either included as a different taxon on the respective taxonomic level or not used for further analysis.

Statistical analysis
All follow-up analyses, as described below, were conducted in R (version 4.0.2, R Core Team, 2020). For significance tests, the significance level was set at p < 0.05.

Metabarcoding of eukaryotic plankton interactions
For alpha diversity, the number of OTUs per season was counted. For principal component analysis (PCA) and beta diversity calculation the OTU table (280 samples with 2790 OTUs) was centred log-ratio (clr) transformed. The PCA was conducted with the "pca" function of the mixOmics package (Rohart et al., 2017). Beta diversity was calculated with the Whittaker index (Whittaker, 1960) using the "betadiver" function of the vegan package (Oksanen et al., 2019). The "hclust" function was used to convert the matrix into a cluster, which was then cut into five branches (at h ¼ 0.8) and visualized. Additionally, clusters were defined at h ¼ 0.5. The clusters at h ¼ 0.5 were tested for significance with a permutational multivariate analysis of variance (PERMANOVA) using the "adonis2" function in the vegan package with 10 000 permutations. This analysis was also used to check for significances of the beta diversity matrix to the different seasons. In order to compare environmental parameters with the defined clusters, we applied the "mantel" function from the ade4 package (Dray and Dufour, 2007;Bougeard and Dray, 2018) with 10 000 permutations. Euclidean distance was used for the dissimilarity matrices of the environmental data and the phases. Missing values in environmental data resulted in deletion of samples before creation of the respective dissimilarity matrices for each analysis (see Supplementary Table S1). We calculated the Shannon and Simpson diversity indices using the "diversitycomp" command of the BiodiversityR package (Kindt and Coe, 2005). Hereby, we used the relative abundances and the seasons of each year and the four seasons combined as grouping variables. Additionally, we visualized the OTU intersection via the "upset" function of the UpSetR package (Conway et al., 2017) for all OTUs per season and for the 200 most abundant OTUs per season. Display of intersections was limited to at least 12 intersections for all OTUs and at least 5 intersections for 200 most abundant OTUs.

Predator-prey interactions and network analysis
As a starting point, we used previously observed predator-prey pairings (Table 1) to check for potential additional information on predator-prey dynamics and grazing impacts. Exemplary successions of main predators and prey were compared, and additional examples of predator-prey pairings were displayed.
Additionally, we inferred planktonic interactions by developing networks for each season that included at least 20 samples. Therefore, summer and autumn 2016 (five and one samples, respectively) and spring 2019 (three samples) were excluded from the analysis. Additionally, we created a network out of the 200 most abundant OTUs of all samples throughout all seasons. The OTU table was converted into relative abundances, prepared as a phyloseq object (McMurdie and Holmes, 2013), and the data were divided into subsets for each season. Due to the uneven number of samples and the high number of OTUs per season, which increased the computational complexity and duration of analysis, we tried to improve the comparison of the single networks by including the same amount of OTUs (200 most abundant each). We then used the method developed by Kurtz et al. (2015) called SParse InversE Covariance estimation for Ecological Association and Statistical Inference (SPIEC-EASI, version 1.1.0) for network construction. In contrast to other methods such as SparCC or Pearson/Spearman, which are based on empirical correlation or co-variance estimations, the SPIEC-EASI method uses the concept of conditional independence. Edges between any two OTUs (network nodes) therefore imply that there exists a relationship between the OTU abundances (association or interaction), which cannot be better explained by creating any other nodes in the network (Kurtz et al., 2015). The pipeline includes a data pre-processing and transformation step (clr transformation). The chosen graphical model was the Meinshausen-Buhlmann's neighbourhood selection, with lambda min ratio set to 1e-2, nlambda set to 20, and 999 repetitions. For visualization of the networks, the results were used to create igraph objects (Csardi and Nepusz, 2006) and plotted with the "plot_network" function of the phyloseq package. The edge width displays edge weights (strength of association). Therefore, positive and negative networks were plotted separately. Node sizes were set proportional to the abundances of the respective OTUs. Additionally, we created a network out of the 200 most abundant OTUs of all samples throughout all seasons using the same parameters as described for the seasonal networks.
These networks may provide further insights into previously observed predator-prey interactions. We tested if new potential Table 1. Potential predator-prey relationships as found by feeding experiments of exemplary taxa.

Predator
Prey or prey preferences References Calanus helgolandicus Chaetoceros pseudocurvisetus, Lauderia borealis, Gymnodinium splendens, Prorocentrum minimum, Skeletonema marinoi, Rhodomonas baltica (Paffenhöfer, 1970(Paffenhöfer, , 1971Schnack, 1979;Lauritano et al., 2011) Centropages hamatus Prefers ciliates over diatoms (Saage et al., 2009 Protoperidinium cf. divergens Copepod nauplii and eggs (Acartia tonsa) ( Jeong, 1994) pairings can be described as they might be visible in the network, for example, due to the formation of subnetworks (small networks consisting of several OTUs, that are not connected to the main network) or the formation of clusters (OTUs of the same taxon that are connected in a main network). Hereby, we assume a positive association between prey and predators is caused by a bottom-up effect, since more prey might lead to more predators. Negative associations are assumed to portray a top-down effect as the higher predator occurrences would cause lower prey abundances.

Successional patterns of different food web components
Especially during spring, we observed a shift from autotrophic Bacillariophyta to Prymnesiophyceae (haptophytes), Trebouxiophyceae, and Ulvophyceae (green algae). Maximum relative sequence abundances of several Bacillariophyta taxa were found in early spring and summer 2016, spring of 2017 and 2018. Highly abundant Bacillariophyta included Rhizosolenia, Thalassiosira, Coscinodiscus, and Ditylum. While the two green macroalgal Ulvophyceae (Ulva and Dilabifilum) were most abundant in spring 2016, spring 2018, and summer 2018, the Prymnesiophyceae (mainly Emiliania and Phaeocystis) constituted more than 10% in spring of 2016 and 2018. The Trebouxiophyceae (Picochlorum) were found both in spring and summer 2018 (for more detailed information on succession of different taxa, see Supplementary Text).
Mixotrophic Dinophyceae shifted from genera such as Gymnodinium and Heterocapsa also occurring during spring to taxa such as Alexandrium in summer to Akashiwo, which mostly occurred in autumn. Other highly abundant genera included Lepidodinium, Tripos, and Prorocentrum, and Chrysochromulina sp. (Prymnesiophyceae).
Heterotroph microzooplankton was mostly represented by Dinophyceae (Gyrodinium sp.), which were most abundant in spring and summer of 2016, 2017, and 2018 and additionally in autumn 2018. Gyrodinium was one of the few genera, which occurred in every sample (Supplementary Table S2 The phylum Metazoa included 20 classes of potential mesoand macrozooplankton taxa, of which 130 genera were identified. The highest relative sequence abundances (i.e. above 10%) were found for different classes of worms (Annelida, Platyhelminthes, Nemertea), for fish (Craniata) as well as for Mollusca, Cnidaria, Ctenophora, and Brachiopoda. Arthropoda was found to be the most abundant class. Arthropoda were present in high relative sequence abundances during all seasons, except for autumn 2016. Brachiopoda (Phoronis) relative sequence abundances were especially high from autumn 2016 onwards until winter 2017/ 2018 and in summer 2018 (Supplementary Table S3, for more detailed information on succession of different taxa see Supplementary Text).

Comparison of planktonic predator occurrence
Most OTUs (out of 2790) represented Opisthokonta and Alveolata (Apicomplexa, Ciliophora, Dinoflagellata, and Perkinsea), 722 and 671 OTUs, respectively (Supplementary Figure 1a). These kingdoms had the highest relative sequence abundances, up to 89.5% and 87.0%, respectively (for details on the succession of other kingdoms see Supplementary Text).
Opisthokonta representatives included consumers such as Metazoa and Choanoflagellida. Decomposers such as Fungi, and Mesomycetozoa were also found ( Table 2). Fungi were highly abundant in summer 2018 with relative sequence abundances of up to 75% (mainly Aspergillus sp. and Cryptococcus sp., see Supplementary Tables S3 and S4), which were co-occurrent with several community shifts. Mesomycetozoa (including, e.g. the parasite Pseudoperkinsus) were most abundant during summers 2016 and 2018 at up to 18% relative abundance.
Out of the most abundant phyla, three groups of planktonic predators were identified: Dinoflagellata, Ciliophora, and Copepoda (Figure 2a  were higher than those of Copepoda and Ciliophora. Copepoda and Dinoflagellata contrasted in relative sequence abundances. If Dinoflagellata relative sequence abundances were high, Copepoda relative sequence abundances were low and vice versa. In two samples (22 February 2018 and 08 March 2018), Ciliophora were more abundant than Dinoflagellata (excluding Syndiniales). In one sample (08 March 2018), Ciliophora had the highest abundance out of the three predator groups and relative sequence abundances of all three predator groups were below 5%. Higher relative sequence abundances of Ciliophora were not only found in occasional samples but throughout several samples of consecutive sampling days, in which Ciliophora had higher relative sequence abundances than Copepoda. This happened for example in May 2016, in May/June 2017, and in March 2018. In general, Ciliophora were less abundant in relation to other predators, although they presented with a high diversity in OTUs (201). In comparison, over 74% of the crustacean OTUs (in total 339) belonged to only two genera (Pseudocalanus and Calanus) and most crustacean OTUs belonged to Copepoda (Supplementary  Table S5). Over 45% of the Dinoflagellata OTUs (in total 442) were parasitic Syndiniales and over 26% remained   connections between the known pairings were observed: several prey preferences are known for the genera Protoperidinium and some distinct species. In our data set three Protoperidinium OTUs were found in low relative sequence abundances (under 1%). Potential prey, such as Skeletonema (Figure 2), had much higher relative sequence abundances and therefore no clear relationships between predator and prey could be assumed. Additionally, even though Thalassiosira and Ditylum (both not displayed) were occasionally co-occurring with Protoperidinium, most peaks were not correlated with the predator. Due to the low relative sequence abundances, no Protoperidinium OTU was part of the network analysis.
Similarly, for copepod predators, such as Paracalanus, Calanus, or Centropages, no clear connection to potential prey could be found (Figure 2). For example, Chaetoceros OTUs were peaking either when Calanus sp. was absent or peaked after the occurrence of predators. However, BLAST alignment did not reveal any Chaetoceros pseudocurvisetus OTU, which was known as prey for Calanus (Table 1). Other potential prey of Calanus and Paracalanus was too low in abundance (Lauderia sp.), only present when the predator was absent (Skeletonema, Rhodomonas) or co-occurring and alternating in peaks (Gymnodinium sp., Prorocentrum sp.) without any distinct patterns. Instead, all investigated relationships revealed high variabilities in occurrences and relative sequence abundances.
We then performed a network analysis to inspect the interconnectivity of the food web during the different seasons and for identification of potential separate networks in the food web. Networks were thus constructed with the 200 most abundant OTUs (based on relative abundance) of each season (see the section Community diversity on comparison of the 200 most abundant OTUs).
We found that the different taxa were highly interconnected with each other in each season (Figures 3 and 4), as OTUs were positively associated with other OTUs of all kinds of taxa at high frequency in one network. A similar structure was also found for the positive association network, which included the 200 most abundant OTUs of all samples ( Figure 5). However, interconnections were varying greatly in strength (thickness of the edges) (Figures 3-5, a list of all associations can be found in Supplementary Table S6). Strength of the association was not associated with abundance of the OTUs, as strong associations were also found between OTUs of high and low relative sequence abundances. Some networks revealed small subnetworks of up to 3 OTUs, which were only positively associated with each other but not to the rest of the main network. These subnetworks did not reveal any separate food web connections; instead subnetworks rather consisted of OTUs of the same taxa, which hinted at these taxa sporadically occurring in high relative sequence abundances. Some OTUs were not found to be positively associated with any other OTU. These included especially OTUs of Opisthokonta, namely of the phyla Metazoa and Fungi, and some Dinoflagellata.
For negative association networks (Supplementary Figure S3) fewer links between different OTUs were found. In spring 2016, the most complex network of negative associations was detected. Overall, most negative associations were found between two single OTUs. Associations occurred between different phyla, but also within single phyla. For example, in winter 2017/2018 six Dinoflagellata OTUs formed a subnetwork, including two parasitic Syndiniales OTUs and four Dinophyceae OTUs.
The network including positive associations throughout all seasons of all samples revealed two subnetworks, which consisted of two OTUs each, as well as three OTUs that were not connected to any other OTU ( Figure 5). One subnetwork consisted of Aspergillus (OTU 52 and 141), the other subnetwork consisted of Temora sp. (OTU 2 and 45). The three separate OTUs all belonged to different Metazoa: OTU 13 (Ctenophora), OTU 29 (Hiatella sp.), and OTU 191 (Metridium sp.). A cluster of eight interconnected Ciliophora OTUs (Choreotrichida and Strombidiida) was found in the main network. This cluster was connected to several other taxa, most connections belonged to OTUs of Chlorophyta, Dinoflagellata, Cryptophyta, and Cercozoa. Additionally, Dinoflagellata OTUs were highly interconnected to each other in the main network and several clusters were formed. For example, several OTUs of Dinophyceae (OTU  21,24,30,36,69,72,85,117,129) and Syndiniales (OTU 54, 100, 163) were highly interconnected, but also connected to further Dinoflagellata OTUs, as well as to other taxa such as Metazoa, Ochrophyta, Chlorophyta, and Ciliophora. Analysis of network connections of OTUs belonging to the observed predator-prey interactions as listed in Table 1 revealed complex structures and a variety of potential predator-prey constellations.
For Paracalanus network analysis revealed positive and negative associations with 10 different phyla: Cercozoa, Chlorophyta, Ciliophora, Dinoflagellata, Fungi, Haptophyta, Metazoa, Ochrophyta, Picozoa, and Stramenopiles_X. Most associations   Calanus OTUs were only present in networks from winter and spring. No negative associations were found. For positive associations, Calanus OTUs were mostly interconnected with other Calanus OTUs. In total, positive associations were found for six phyla: Ciliophora, Cryptophyta, Dinoflagellata, unclassified Eukaryota, Metazoa, and Ochrophyta. Potential prey as observed before (Table 1) did not show any connections. The only association to a diatom OTU (OTU 307) was found in the winter 2017/ 2018 network, displaying the strongest connection besides the interconnections of different Calanus OTUs. In winter 2018/ 2019, OTU 123 that was identified as Cryptomonadales was connected to two Calanus OTUs (OTU 333 and OTU 593). In spring 2018, a connection to OTU 306 (unclassified Gymnodiniales) and an even stronger connection to a Chrysophyceae OTU (OTU 844) were found.
Centropages was associated with 7 different phyla: Cercozoa, Choanoflagellida, Ciliophora, Dinoflagellata, Metazoa, Ochrophyta, and Stramenopiles_X. Most associations were positive, only one negative association to a ciliate (OTU 248) was found in winter 2017/2018. As indicated by Table 1, Centropages seems to prefer ciliates over diatoms, as more connections to different ciliate OTUs were found compared to connections to diatoms. In terms of strength of association, connections to ciliates were stronger in 2017 compared to 2018. Positive associations of ciliates were found in summer 2017, autumn 2017, and autumn 2018, for OTU 497, OTU 130, and OTU 590, respectively. A positive association to diatom OTU 487 was found in spring 2018, a weaker association to diatom OTU 84 in summer 2018.

Community diversity
Community composition significantly changed from one season to the next during all three seasonal cycles, while the communities of the individual seasons observed in the different years were highly similar. This is reflected in the PCA plot ( Figure 6) showing that samples from the same season rather than samples from the same year cluster together. Samples collected in spring 2017 showed the highest inter-sample variability compared to other years, whereas spring samples were, in general, more similar to each other than autumn samples. The two outliers of the PCA plot from the 01 February 2018 (winter 2017/2018) and 18 December 2018 (winter 2018/2019) can be explained by the heterogeneity of library sizes, since these two samples consisted of a bigger library than all other samples.
In most years, the diversity of the plankton communities was higher in autumn and winter than in spring and summer, while summer displayed the lowest diversity. This is reflected by both, OTU numbers (Figure 7, Supplementary Figure S4) and diversity indices (Tables 3 and 4), whereas Shannon and Simpson diversity indices were not differing much in general. Sampling intensity for the different seasons of the different years was variable due to logistic constraints. For most of the seasons, more samples collected for a respective season did not result in a higher diversity (sample size above 20). For example, even though the number of samples taken in spring 2016 was nearly twice as much as in spring 2017 and spring 2018, 49 samples compared to 25 samples each, the number of OTUs was similar, 2098, 2092, and 2016 OTUs, respectively. Considering Shannon and Simpson indices, species composition of autumn and winter samples was most diverse, while spring diversity was lower and summer communities had the lowest diversity (Table 3). However, the subset of three seasons with the fewest number of samples, also had the least number of OTUs as well as the lowest Shannon and Simpson diversity, whereas both indices resulted in high values in general (Table 4)

Metabarcoding of eukaryotic plankton interactions
during spring more OTUs of the class Cercozoa belonged to the 200 most abundant taxa than during other seasons. In total, our data set for analysis consisted of 2790 OTUs. Between 2104 and 2313 OTUs were found during winter. Alpha diversity in autumn was between 834 and 2156 OTUs. Out of all OTUs, 295 OTUs (10.6%) were present throughout all seasons (Figure 7). Furthermore, 168 OTUs were shared by all seasons, except for autumn 2016, which were sampled just once. Only 13 OTUs that belonged to the most abundant OTUs were present in all seasons (Supplementary Figure S4). A higher number of unique OTUs (Figure 7) and the highest proportion of unique Figure  S4) indicated that the community of summer 2018 was different from other years and seasons. The beta diversity analysis further indicated that spring and summer 2018 were different from other years, as the samples from these seasons resulted in several significant small clusters (Supplementary Figures S5 and S6, more details in Supplementary Text, R 2 ¼0.69727, F ¼ 16.597, p < 0.0001). PERMANOVA also confirmed that the communities occurring during each season (13 seasons, R 2 ¼0.47509, F ¼ 20.138 p < 0.0001) and the OTUs of the four seasons across all years (R 2 ¼0.31119, F ¼ 41.564, p < 0.0001) were significantly different. While this difference was also caused by the different environmental conditions, the Mantel test revealed a significant but mostly weak correlation of the metabarcoding data set with several environmental factors: temperature, nitrate, sunshine duration, salinity, and Secchi depth (see Supplementary Table S1).

Discussion
Our 3-year metabarcoding study revealed a highly variable system, in which blooms of single prey taxa are only occurring occasionally and without any distinct pattern, whereas potential predators are found in high relative sequence abundances throughout. While the overall abundance of Bacillariophyta was highest in spring followed by summer and autumn, certain genera did not bloom during every year. Throughout the years, similar findings on changes in abundances of diatoms, but also of shifts in occurrence were recorded (Scharfe and Wiltshire, 2019).
In the following paragraphs, we, therefore, want to discuss the following important results: We found a very strong inter-annual variation in the algae, but not in the grazers and existing predator-prey relationships could not be found in the metabarcoding data. Instead, our networks show very strong connections with many nodes, indicating that the webs are probably very robust, and the predators seem to be able to shift without any problems from one prey item to another.

Community diversity
Using water samples, we detected a high diversity of taxa, ranging from potential meroplankton, such as fish and other Metazoa, to   holoplanktonic mesoplankton down to the smallest picoplankton and parasites. Parasitic taxa, which mostly are represented by marine parasitoids, hereby can include free-living stages, but also parasitoids currently infecting other plankton (Käse et al., 2021). Intensive sampling revealed a high diversity and likely increased the probability of catching rarer taxa, as confirmed by both, Shannon and Simpson, diversity indices. However, a maximum diversity was reached as more sampling did not result in a higher diversity, e.g. when comparing sampling efforts during different spring seasons.
Furthermore, despite spring blooms being frequently monitored, spring revealed less OTUs than other seasons, especially compared to winter and resulted in lower diversity indices. Instead, autumn was the most diverse season. Comparison of OTUs per season revealed a steady background community consisting of various taxa (295 OTUs), which were sampled during every season, despite the strong seasonality at the sampling site. The main taxa on phylum level were also, with few exceptions, present nearly all the time. A diverse background community of protists was reported in spring 2016 , which therefore can now be extrapolated to other seasons and includes also bigger sized zooplankton.
However, the presence of a seemingly stable background community, which has now been shown in two studies, does not mean considerable fluctuations such as blooms of unusual species are not possible. This is exemplified by the year 2018. This year was unique in terms of community composition, with summer samples, in particular, differing markedly as seen in the amount of unique OTUs. The occurrence of fungi, Picochlorum and Dilabifilum hints at a community of the algae and lichen-forming fungi which was previously observed for different Ulvophyceae and lichen-forming fungi Thüs et al., 2011). In addition, benthic taxa, such as worms, were occasionally found in high relative sequence abundances. There are several potential explanations for this. The most parsimonious is that due to the relatively shallow sampling site, combined with the strong tidal influence, and the occasional storm, material, and organisms were suspended into the water column and sampled in our water samples. However, it is also possible that these species were sampled in their (mero-)planktonic state instead of the adult state, which is indistinguishable by metabarcoding (Bucklin et al., 2016).

Predator-prey interactions and network analysis
We did not observe any strong support of the predator-prey pairs that we know, nor did we find any close connections to others. One reason for this might be our observation that even though the prey communities rapidly change and do not re-occur with the same species from 1 year to the next, this is completely different for predators. Hence, a strong link between individual taxa cannot be expected.
Dinoflagellates made the highest contribution to the community, most likely, playing a key role not only as a predator but also as potential prey for bigger sized taxa. A general bias of our primers in favour of dinoflagellates is unlikely, even though dinoflagellates have a higher copy number, as Sprong et al. (2020) showed that more coastal stations were not dinoflagellate dominated using the same primers compared to our sampling area. Similar results were found in the Estuary of Bilbao, where high relative sequence abundances of copepods were found in larger size fractions, and no dominance of dinoflagellates even in small size fractions was seen (Abad et al., 2017). Therefore, our results were most likely a reflection of the unique ecology at our sampling point and not caused by a bias in the molecular method.
The high relative sequence abundances of heterotrophic dinoflagellates are supported by a previous study, using traditional microscopy, which also detected high contribution of heterotrophic dinoflagellates in biomass (Löder et al., 2012). Our observation of ciliates peaking during spring is also supported by previous studies (Löder et al., 2011(Löder et al., , 2012, which found ciliates to be an important but highly selective grazer in spring. Besides the possibility of methodological constraints regarding ciliate detection, this specificity could explain the low relative sequence abundances in our assemblage during certain years. The variability of diatom occurrences as prey might be reflected by the grazer relative sequence abundances as well. However, a potential cluster of ciliates was found by network analysis, and a variety of taxa were associated with this cluster, indicating several potential predatorprey relationships. Furthermore, it needs to be noted that bacteria, which were not part of this study, are known as an additional important prey option for ciliates (Albright et al., 1987;Sherr and Sherr, 1987).
The copepods Paracalanus and Centropages, which are able to feed on ciliates (Suzuki et al., 1999;Saage et al., 2009), were associated with ciliates in the network analysis. However, in general, connections were weak, and stronger connections to other taxa were found. In contrast to ciliates, crustacean OTUs peaked every year no matter which potential food was present. While copepods can feed selectively on microzooplankton (Nejstgaard et al., 2001;Löder et al., 2011), they are often omnivorous, feeding on a diverse range of organisms such as diatoms, dinoflagellates, and other zooplankton including their own eggs and larvae (e.g. Calbet and Saiz, 2005;Boersma et al., 2014;Yeh et al., 2020).
The difficulty of identifying predator-prey dynamics could be explained by potential coexistence within different plankton communities due to different feeding and motility traits as it has been shown, e.g. for the western English Channel, (Kenitz et al., 2017). The authors also linked seasonal succession in the community, besides the interannual variation in dominant species, to a succession of activity traits. Furthermore, Kenitz et al. (2017) suggested that strong turbulence benefit passive feeding zooplankton, and leads to an enhanced grazing pressure on motile prey, which would benefit the growth of non-motile prey. In our data set, Oithona, as a passive feeder, occurred in high relative sequence abundances in summer but was also abundant during other seasons and we observed additional blooms of non-motile prey during summer months, which might have benefitted from decreased grazing pressure.
Similar to the differences in copy numbers for dinoflagellates, a large proportion of DNA of multi-celled zooplankton such as copepods or fish might indicate a bias of these taxa in the metabarcoding assemblage. However, in our approach a dominance of these taxa was not evident, as generally, dinoflagellates were most abundant, nor did high relative sequence abundances prevent the formation of seasonal patterns of small-sized plankton. This further indicates that the data on these taxa was mostly based on DNA, which was not retrieved through extraction of whole organisms or their extremities, but rather from environmental DNA. Due to continuously high relative sequence abundances and diversity of the predator and prey as well as the complexity of the food web, a distinct grazing impact on single taxa or distinct links between potential predator-prey pairings could not be extracted from our data set. While it is known, that especially zooplankton biomass responds in much longer time scales , e.g. due to the complexity of the metazoan life stages, we pose that co-occurrence might already be hinting at a potential relationship. For example, Calanus might prey on Chaetoceros or peaks of predators might indicate that prey is eaten up and therefore no longer detected. It could also be the case that peaks of predators are caused by feeding on other prey instead, which makes it difficult to define distinct predatorprey relationships. However, other investigated potential distinct predator-prey relationships could also not be observed.
In general, network analysis of the 200 most abundant OTUs revealed Metazoa OTUs, whose respective organisms would be bigger in size than the rest of the sampled plankton community, were not as tightly connected to the rest of the network. This might be the case, especially if potential consumers are highly abundant on rare occasions. In general, these networks could be a potential tool for detection of specific relationships. However, the networks were tightly interconnected and only few OTUs did not connect to the main seasonal networks. Links between edges of the same taxonomic group have been observed before (Faust et al., 2012;Kurtz et al., 2015) and are commonly described as "assortativity." These associations could also be detected because occasionally two OTUs represent the same species but could not be merged to a single OTU because of potential errors in the sequence or other biases. Additionally, the risk of spurious coincidences might be increased in the network analysis and associations might be depicted by chance. This is why interpretation needs to be careful and additional analyses are necessary to verify potential associations. As each season also comprises several communities as depicted by the clusters in the beta diversity, the networks might also include associations in between these different communities. Especially since the PCA analysis also indicated that few samples might rather belong to a previous or follow-up season, more or other associations might have been found without the focus on the seasonality. Based on these complex networks, clear dynamics cannot be identified easily and interactions between food web components seem to occur on much bigger scales. The tight links in between various components of the food web emphasize that they are co-occurring throughout the different seasons and indicate a high variability in food options for potential consumers.
Previous investigations of copepod faecal pellets (Turner, 1984) and metabarcoding of the gut content (Yeh et al., 2020) revealed that copepods ingest a wide variety of food items. Besides this high variability in ingested food, combining the known predator-prey pairings of previous grazing experiments also demonstrates this high variability. The selective feeding on certain taxa has mostly been observed in grazing experiments, where potential prey is limited to certain taxa and provided constantly, whereas the complex hydrography at Helgoland might disturb the actual formation of a system sufficiently stable for allowing specific predator-prey relationships. Therefore, it is possible that predators are more flexible and less selective in natural environments than in experiments, but they are also provided with a higher choice and variability of prey options. Alternatively, we are not able to distinguish existing predator-prey relationships as the high variability in the system conceals explicit dynamics.
The fact that we found a background community in addition to the very variable occurrence of taxa, which might include mostly opportunistic species, hints at a rather flexible but nevertheless stable and healthy ecosystem. While shifts of species occurrence due to environmental change were observed at Helgoland (Scharfe and Wiltshire, 2019), stability of the ecosystem might be high enough due to the natural fluctuation and high adaptability of the system. The tight links in the seasonal networks indicate furthermore, that the robustness of the food web to species loss is potentially quite high (Dunne et al., 2002;Estrada, 2007), which is also supported by the random occurrence of taxa throughout the years.

Conclusion
Metabarcoding of water samples is suitable for capturing taxa of the whole community and for obtaining information on plankton succession in relation to time and environmental conditions, without exclusion of large-sized taxa from the analysis. Therefore, new technologies, like next-generation sequencing, may be used in addition to traditional methodologies on a long-term basis. Comparability and practicability of combining these different methods still need to be tested in future studies. The system is characterized by a high variability of potential prey or predators, which are not necessarily co-occurring or displaying typical patterns. Predator-prey relationships in the planktonic community are diverse and plentiful and specialist relationships are rather uncommon. Instead, generalists seem to be the norm, which makes it difficult to extract distinct predator-prey dynamics in the field. This offers an enormous potential of relationships in the plankton community that might be verified by traditional laboratory experiments. Furthermore, it remains under question how the resilience of the North Sea is influenced by the high variability.

Supplementary data
Supplementary material is available at the ICESJMS online version of the manuscript.

Data availability
Sequence data for this study have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB37135, using the data brokerage service of the German Federation for Biological Data (GFBio, Diepenbroek et al., 2014), in compliance with the Minimal Information about any (X) Sequence (MIxS) standard (Yilmaz et al., 2011). Details of our pipeline are available on GitHub (https://github.com/PyoneerO/ qzip). Additionally, the full OTU table (280 samples with 59.284 OTUs) was archived in PANGAEA (DOI: 10.1594/ PANGAEA.921026). All other relevant data are within the manuscript and its Supplementary material files.