Meroplankton Diversity, Seasonality and Life-History Traits Across the Barents Sea Polar Front Revealed by High-Throughput DNA Barcoding

In many species of marine benthic invertebrates, a planktonic larval phase plays a critical role in dispersal. Very little is known about the larval biology of most species, however, in part because species identification has historically been hindered by the microscopic size and morphological similarity among related taxa. This study aimed to determine the taxonomic composition and seasonal distribution of meroplankton in the Barents Sea, across the Polar Front. We collected meroplankton during five time points seasonally and used high-throughput DNA barcoding of individual larvae to obtain species-level information on larval seasonality. We found that meroplankton was highly diverse (72 taxa from eight phyla) and present in the Barents Sea year-round with a peak in abundance in August and November, defying the conventional wisdom that peak abundance would coincide with the spring phytoplankton bloom. Ophiuroids, bivalves, and polychaetes dominated larval abundance while gastropods and polychaetes accounted for the bulk of the taxon diversity. Community structure varied seasonally and total abundance was generally higher south of the Polar Front while taxon richness was overall greater to the north. Of the species identified, most were known inhabitants of the Barents Sea. However, the nemertean Cephalothrix iwatai and the brittle star Ophiocten gracilis were abundant in the meroplankton despite never having been previously recorded in the northern Barents Sea. The new knowledge on seasonal patterns of individual meroplanktonic species has implications for understanding environment-biotic interactions in a changing Arctic and provides a framework for early detection of potential newcomers to the system.

In many species of marine benthic invertebrates, a planktonic larval phase plays a critical role in dispersal. Very little is known about the larval biology of most species, however, in part because species identification has historically been hindered by the microscopic size and morphological similarity among related taxa. This study aimed to determine the taxonomic composition and seasonal distribution of meroplankton in the Barents Sea, across the Polar Front. We collected meroplankton during five time points seasonally and used high-throughput DNA barcoding of individual larvae to obtain species-level information on larval seasonality. We found that meroplankton was highly diverse (72 taxa from eight phyla) and present in the Barents Sea yearround with a peak in abundance in August and November, defying the conventional wisdom that peak abundance would coincide with the spring phytoplankton bloom. Ophiuroids, bivalves, and polychaetes dominated larval abundance while gastropods and polychaetes accounted for the bulk of the taxon diversity. Community structure varied seasonally and total abundance was generally higher south of the Polar Front while taxon richness was overall greater to the north. Of the species identified, most were known inhabitants of the Barents Sea. However, the nemertean Cephalothrix iwatai and the brittle star Ophiocten gracilis were abundant in the meroplankton despite never having been previously recorded in the northern Barents Sea. The new knowledge on seasonal patterns of individual meroplanktonic species has implications for understanding environment-biotic interactions in a changing Arctic and provides a framework for early detection of potential newcomers to the system.

INTRODUCTION
In the Arctic, and around the world's oceans, benthic invertebrates play important roles in carbon cycling and remineralization of nutrients (Renaud et al., 2007a), as prey for fish (Eriksen et al., 2020), birds (Merkel et al., 2007), and mammals (Dehn et al., 2007), as well as supporting important fisheries (e.g., the Northern shrimp Pandalus borealis, Garcia, 2007), and subsistence harvesting (Rapinski et al., 2018). The Barents Sea, located in the Atlantic gateway to the Arctic, is home to over 3,000 benthic invertebrate taxa, making it one of the most diverse regions of the Arctic (Piepenburg et al., 2011;Renaud et al., 2015).
Most studies on the benthic invertebrates of the Barents Sea have focused on the adult stage (Carroll et al., 2008;Cochrane et al., 2009;Jørgensen et al., 2015;Zakharov et al., 2020). For a large proportion of benthic invertebrates, however, the larval phase constitutes the main vector for dispersal, having important implications for population connectivity (Cowen and Sponaugle, 2009), species conservation and management of marine protected areas (Krueck et al., 2017), dispersal of invasive species and understanding of the impacts of environmental change on benthic communities (Levin, 2006;Renaud et al., 2015). Compared to adults, planktonic stages are relatively shortlived, with a presence in the water column from hours to months (Shanks, 2009), though this becomes longer in colder temperatures, sometimes more than doubling in duration with a drop of a few degrees (O'Connor et al., 2007). Therefore, repeated sampling over the course of a year is required to capture as much of the diversity as possible. Such seasonal sampling of Arctic meroplankton has mostly been done in fjordic environments to date (Kuklinski et al., 2013;Stübner et al., 2016;Brandner et al., 2017;Michelsen et al., 2017). There, most meroplanktonic groups peak in abundance in the spring/early summer, around the time of the phytoplankton bloom, presumably to maximize food availability for the larvae. The few studies that focus on Arctic offshore locations only sampled during a short time window, together covering only the period from May to September (Clough et al., 1997;Schlüter and Rachor, 2001;Fetzer and Arntz, 2008). A peak in larval abundance coinciding with a peak in local primary production, as occurs in Arctic and Antarctic fjords and coasts (Bowden et al., 2009;Arendt et al., 2013;Michelsen et al., 2017;Presta et al., 2019) as well as in lower latitude regions (Highfield et al., 2010), is often assumed but has not yet been corroborated with seasonally resolved sampling on Arctic shelves.
The spring phytoplankton bloom is broadly found to occur in May in the southern Barents Sea, when sufficient light and stratification of the water column favor bloom development. In ice-covered waters of the northern Barents Sea, however, the phenology of the phytoplankton bloom is more variable, occurring anytime from May to July depending on timing of sea ice retreat (Dalpadado et al., 2020;Dong et al., 2020). In these seasonally ice-covered waters, an ice algal bloom as well as an under-ice phytoplankton bloom can contribute primary production to the system prior to ice melt, extending the duration of the productive period (Syvertsen, 1991;Leu et al., 2015;Ardyna et al., 2020). Even with these additional sources of production in ice-covered waters, however, total primary production in the northern Barents Sea is substantially lower than in the southern Barents Sea (Slagstad et al., 2011;Hunt et al., 2013).
The Barents Sea is an Arctic inflow shelf which receives large inputs of Atlantic Water (Carmack and Wassmann, 2006;Jakobsen et al., 2011). As this Atlantic Water flows north and eastward, it eventually meets colder and fresher Arctic Water flowing south (Oziel et al., 2016; Figure 1). The boundary area between these two major water masses is termed the Polar Front (Oziel et al., 2016). The front is relatively narrow and its location constrained largely by topography, remaining relatively constant through time. East of approximately 32 • E the Polar Front splits into a southern and a northern front, the locations of which are typically more variable than the front west of 32 • E (Oziel et al., 2016), though Barton et al. (2018) noted that the location of the northern front is also well constrained by topography. As Atlantic Water cools and mixes with other water masses in the Barents Sea, it becomes denser and sinks, forming Barents Sea Water (Oziel et al., 2016) which can then penetrate under Arctic Water north of the Polar Front (Lind et al., 2012). Atlantic Water entering the Barents Sea from the south can supply high abundances of biota, including zooplankton, onto the Arctic shelf (Wassmann et al., 2015;Hunt et al., 2016), potentially including planktonic life stages of more southern benthic species. Arctic inflow shelves such as the Barents Sea are thus likely to increasingly receive boreal species extending their range into the Arctic as the climate continues to warm (Renaud et al., 2015). At this time, we do not know which species have the capacity to reach the Arctic during their larval stage.
Benthic invertebrate early life stages can be very difficult, and sometimes even impossible, to identify to species or genus using morphological characteristics alone. Consequently, most zooplankton studies either group taxa into an all-encompassing "meroplankton" category or identify them to broad taxonomic resolution such as phylum or class levels (e.g., Hirche et al., 2011;Gluchowska et al., 2016). Even in studies that focus specifically on meroplankton, species-level identification is scarce due to the common lack of distinguishing features across closely-related taxa (e.g., Bowden et al., 2009;Highfield et al., 2010). The use of DNA barcoding to identify meroplankton, while promising, initially achieved relatively low success rates (Heimeier et al., 2010). However, recent advances in methods specifically designed to barcode meroplankton have enabled the identification of hundreds of individuals, and dozens of taxa, with high success rate (Ershova et al., 2019).
Such species-level identification opens up possibilities to obtain detailed community composition and its variability in time and space. In addition, species-level identification can contribute important information on larval biological traits such as timing of larval release as well as larval duration for a large number of species simultaneously. Historically, much of the information on larval duration of individual species relied on laboratory rearing of larvae (Shanks, 2009). As such, the information was limited to a subset of species, each experiment studying a single or a few species at once, and could have potentially been biased by laboratory conditions (Shanks, 2009). Finally, species-level identification also has the potential to detect non-native species as they drift into the Arctic through water currents, perhaps serving as an early warning system for harmful alien species, or species of potential ecosystem impact.
The objective of this study was to determine the taxonomic composition and seasonality of the meroplankton community in the Barents Sea, north and south of the Polar Front. By using a high-throughput DNA barcoding method to identify early life stages, we aimed to achieve high taxonomic resolution yielding species-specific information on seasonality, larval duration, and FIGURE 1 | Study area with sampling locations. Map was produced using package ggOceanMaps in R (Vihtakari, 2020). Arrows show Atlantic (AW, red) and Arctic (ArW, blue) water currents while the whole dark line represents the approximate location of the oceanographic Polar Front (PF) where it is stable, west of 32 • E, and dashed lines represent the more variable southern and northern fronts to the east (based on Loeng, 1991;Oziel et al., 2016). association with water masses, as well as identifying greater taxon richness than previously realized, including identification of potential newcomers. We hypothesized that seasonality differs across taxa but that most species peak in abundance around the spring phytoplankton bloom and remain in the water column for up to a few months. We also predicted that the meroplankton community differs in composition between north and south of the Polar Front, in association with different water masses, and that the larval assemblage represents a mix of larvae from local Barents Sea species as well as some more southern species advected into the Barents Sea. Implications of these findings in the context of invasibility of the Arctic inflow shelves are discussed.

Field Sampling
Meroplankton samples were collected at stations north and south of the Polar Front during a series of seasonal cruises in November 2017 and in January, April, June and August 2018 (Table 1 and Figure 1). Due to cruise logistics, no samples were collected south of the Polar Front in November 2017. In addition, extensive sea ice blocked northward travel in April, so sampling north of the Polar Front occurred at slightly lower latitude in April compared to the other time points. At each station and sampling month, a CTD cast (Seabird SBE 911 or SBE 911+) preceded zooplankton sampling to determine the depth distribution of water masses based on the vertical profiles of temperature, salinity and density (Table 1). Temperature, salinity, and fluorescence measurements (as a proxy for chlorophyll a concentration) data were obtained from the CTD casts for further analysis. Zooplankton samples were collected from each water layer individually using a 64-µm closing WP2 (0.255 m 2 mouth opening, Hydro-Bios) or Multinet (0.25 m 2 mouth opening, Hydro-Bios) towed vertically at 0.5 m/s, retrieving two to four samples per station and covering the entire water column from ∼10 m above the seafloor to the surface. In April, when phytoplankton was highly abundant, 180 µm mesh was used instead of 64 µm to prevent clogging. The 180 µm mesh was also used to sample the surface layer in June at the northern location after the upper portion of the 64 µm mesh ripped open during sampling of lower layers, leaving only the outer 300 µm mesh toward the net opening (see Table 1 for further details). Zooplankton samples were concentrated on a 64 µm sieve and

Sorting and Quantification
In the laboratory, a quantitative subsample from each zooplankton sample was sorted under a Leica M205C dissecting microscope at 20-100 × magnification to quantify and isolate individual meroplankton specimens for barcoding. Subsamples of 0.8-23.5% of total sample volume ( Table 1), depending on meroplankton density, were collected using a 5-mL pipette with the tip cut off to prevent clogging from large zooplankton. All meroplanktonic individuals in the subsample were visually classified into morphotypes largely corresponding to broad taxonomic levels (e.g., Bivalvia, Nudibranchia, Ophiuroidea, etc.) and counted. For each morphotype, the number of individuals counted in the subsample was multiplied by the subsampling factor to estimate the numbers in the whole sample. Abundance values (individuals/m 3 ) of each morphotype in each sample were then calculated by dividing the number of larvae by the volume of seawater filtered for each sample. For each sample, up to ∼30 individuals from each morphotype were photographed (for reference and for measurements of body size) using a Leica MC170 HD microscope camera, then soaked in MilliQ water prior to DNA extraction.

High-Throughput DNA Barcoding
Each photographed individual was then transferred individually into 20 µL Alkaline Lysis Reagent (25 mM NaOH, 0.2 mM EDTA) using sterilized tools, and heated at 95 • C for 30 min to extract the DNA. Bivalve and gastropod larvae were crushed with the flat blade of a micro-scalpel prior to extraction. The DNA extracts were then pH-neutralized by adding 20 µL Neutralization Buffer (40 mM Tris-HCl) and kept at 4 • C until amplification. This extraction method (HotShot extraction, Truett et al., 2000;Meissner et al., 2013) is ideal for very small organisms like meroplankton because it minimizes DNA loss by eliminating all cleaning and transfer steps. We used a high-throughput multiplexing approach (after Ershova et al., 2019) to DNA-barcode a large number of larvae (Shokralla et al., 2015). By using a double-tagging strategy, we were able to sequence ∼1,000 individuals simultaneously. First, a ∼313 base pair fragment of the mitochondrial cytochrome C oxidase I (COI) gene was amplified by PCR using individually-tagged Leray-XT primers, a highly degenerate primer set (Geller et al., 2013;Wangensteen et al., 2018). The availability of reference sequences as well as its ability to assign taxonomy to species level (Hebert et al., 2003;Andújar et al., 2018) make this marker ideal for meroplankton identification. Each PCR reaction consisted of 10 µL polymerase (AccuStart II PCR ToughMix or AmpliTaq Gold), 0.16 µL bovine serum albumin 20 µg/µL, 5.84 µL nuclease-free water, 1 µL individually-tagged forward primer (5 µM, mlCOIintF-XT 5 -GGWACWRGWTGRACWITITAYCCYCC-3 ), 1 µL individually-tagged reverse primer (5 µM, jgHCO2198 5 -TAIACYTCIGGRTGICCRAARAAYCA-3 ) as well as 2 µL undiluted DNA template for a total reaction volume of 20 µL.
The PCR protocol consisted of a denaturation step for 10 min at 95 • C followed by 35 cycles of 94 • C for 1 min, 45 • C for 1 min and 72 • C for 1 min, and a final extension of 5 min at 72 • C (Wangensteen et al., 2018). The same tag was used on the forward and reverse primers to enable detection of chimeras. A total of 96 different primer tags were available, so sets of 96 individually-tagged PCR amplicons were pooled into libraries. Libraries were cleaned (fragments below 70 base pairs removed) and concentrated using MinElute columns and final DNA concentration measured using a Qubit fluorometer with broad-range dsDNA BR Assay Kit (Qubit). Each library was then prepared for sequencing using the NEXTflex PCR-free DNA sequencing kit (BIOO Scientific), according to manufacturer instructions, omitting the second bead cleaning round in Step B to improve yield. This library preparation protocol includes ligation of Illumina adapters as well as a 6-base library tag (BIOO Scientific). Each pool of 96 samples received a different library tag so that DNA from each larva was recognizable through its unique combination of the primer tag and the library tag. Prior to sequencing, each library was quantified by qPCR using the NEBNext Library Quant Kit (New England Biolabs). Libraries were sequenced on an Illumina MiSeq platform with v2 or v3 2 × 250 bp kits and spiked with 1% PhiX used as an internal control to calculate error rates.

Bioinformatics
All sequences obtained from the MiSeq runs were processed in the OBITools software suite (Boyer et al., 2016) as in Ershova et al. (2019). The forward and reverse reads were aligned and annotated for reads with alignment quality scores > 40 using the function illuminapairedend. Contigs were assigned to the corresponding larvae (demultiplexed) using ngsfilter, which also removed primer sequences. Only sequences with 300-320 base pairs and no ambiguous bases were retained for clustering in Swarm v2 (Mahé et al., 2015) with local clustering threshold d = 13. Singletons were discarded. Initial taxonomic assignation was performed with Ecotag (Boyer et al., 2016) against a local reference database (Wangensteen et al., 2018) and then checked manually in the Barcode of Life database (BOLD 1 ) and NCBI database BLAST 2 for definitive assignment. Due to the nature of high-throughput sequencing, most larval samples contained reads from several taxa (possibly food items, contamination, etc.). The sequence with the most abundant reads within each sample was assumed to correspond to the larva itself. All samples with either less than 10 reads total or with fewer than 1,000 reads and, at the same time, less than 75% of reads assigned to the most abundant sequence were considered unreliable and were not included in further analysis. We then compared the taxonomic assignments from DNA barcoding to the initial morphological identification. Individuals, for which the two did not match at even a coarse taxonomic level, as well as non-meroplanktonic taxa, were omitted from subsequent analysis. The remaining sequences were considered successful. A sequence match of >98% to the reference database was considered an "exact" match (Leray et al., 2016). For the larvae with <98% sequence match, the taxon name of the closest match was retained, but the % match was listed alongside for clarity. Accepted taxonomic names and classification were obtained from the World Register of Marine Species (WoRMS 3 ).

Larval Body Size Measurements
For species that were detected in the samples at multiple time points, body size was measured to help differentiate between species with continuous reproduction (average body size relatively constant throughout the year) from species with more discrete reproductive timing but long larval duration (average body size increases through time). For the purpose of this analysis, temperature effects were not considered. Body size measurements were made using the photographs of each individual obtained prior to DNA extraction. Measurements were conducted in ImageJ software (Abramoff et al., 2004) with the Straight Line tool (Segmented Line when needed), using the scale bar on each photograph to calibrate the measurements. In echinoderm pluteus larvae, the length of the postero-dorsal arm was used while in bivalve veliger larvae, the anterior-posterior length was measured (Shanks, 2001).

Data Analysis
Temperature and salinity data from the station CTD casts were manually de-spiked and binned by 1 m (or 1 db). Salinity data were calibrated to Autosal salinities when available (June-July and August 2018). Each sample from a given water layer was assigned to a water mass based on its mean salinity (S) and temperature (T) values: Atlantic Water (AW: S > 34.8, T > 3 • C), Arctic Water (ArW: S < 34.7, T < 0 • C) and Barents Sea Water (BSW: S > 34.8, T ≤ 2 • C, Oziel et al., 2016). Where mean temperature and salinity values fell between the defined cores of two water masses, the water was considered to be a mixing product of the two. When the surface layer did not fit in the standard water mass definitions but was much fresher than underlying layers, it was called Surface Water (SW).
Density values for each species detected by DNA barcoding were obtained by multiplying the number of each morphotype in each sample by the percentage of the barcoded larvae corresponding to each sequence-based identification. For example, if a sample contained 50 ophioplutei/m 3 and from that sample 10% of all barcoded ophioplutei matched with Ophiura sarsii, then we estimated the density of O. sarsii in that sample to be 5 individuals/m 3 . We then multiplied the density in each sample by the depth range of that same sample and added all depth layers within one location to obtain a depth-integrated abundance value expressed in individuals/m 2 . Larvae that were not successfully identified through DNA barcoding were not included in subsequent analysis so listed abundances represent a conservative estimate of total abundances. Meroplankton diversity and community analyses were conducted with package vegan in RStudio, v. 1.4.1106 (Oksanen et al., 2019) and visualized using package ggplot2 (Wickham, 2016). A species accumulation curve (function specaccum) was used to determine 3 www.marinespecies.org/ whether our sampling covered the full diversity of the region. A heatmap of species abundances was created to visualize seasonal patterns. Meroplanktonic taxa in the heatmap were clustered according to the seasonal and geographical distribution using the UPGMA (unweighted pair group method with arithmetic mean) clustering method. To test the difference in meroplankton communities across months and locations, we did separate analyses of "size" (total abundance) and "shape" (community composition, Greenacre, 2017). The difference in total abundance across locations and months was tested using Kruskal-Wallis tests performed separately for location and month. The unequal variance across locations and months prevented the use of an ANOVA. The difference in taxon richness across location and month, however, was tested using a two-way ANOVA without an interaction term. Spatial and seasonal patterns in community structure were explored using Canonical Correspondence Analysis (CCA) of the fourth-root transformed taxa abundances. Explanatory variables included water mass type (Atlantic, Arctic, Barents Sea and Surface Waters; Table 1), location (north vs. south of the Polar Front), sampling month (November, January, April, June, and August), as well as average water column temperature, salinity, and fluorescence (as a proxy for chlorophyll a concentrations). The best subset of constraining variables was selected by stepwise selection using the ordistep function. The significance of the overall model and individual terms was calculated using permutation tests (anova.cca function) at a significance level of p < 0.05. The correlations of the individual species abundances to the resulting ordination were obtained using the function envfit, and only species that were found to be significantly driving the ordination (p < 0.05) are shown on the biplots.

Oceanography
Water mass characteristics contrasted sharply between the north and south of the Polar Front during all seasons (Figure 2 and Table 1). North of the Polar Front, deep waters consisted of the relatively cold and saline Barents Sea Water throughout the year while shallower layers consisted mostly of Arctic Water. The characteristics of the surface layer fluctuated seasonally. A colder and fresher water layer overlaid the Arctic Water in November, but was not present in January. Surface temperatures began to increase again in June and reached upward of 4 • C by August. While sampling north of the Polar Front in April occurred further south than in the other months (Figure 1), the water mass characteristics showed a clear similarity to the more northern sampling location of other months, remaining well within characteristics of Arctic Water at the surface (Figure 2). The deep waters south of the Polar Front also consisted mostly of Barents Sea Water, except in January where they were typically 1-2 • C warmer. Upper layers in the south consisted entirely of Atlantic Water, except in April when the whole water column was well mixed with characteristics of the Barents Sea Water mass. Together, these observations suggest that our two sampling locations were positioned on either side of the FIGURE 2 | Temperature and salinity diagram of water masses in the Barents Sea where meroplankton was collected. The colored rectangles represent the major water masses encountered during this study: Atlantic Water (red), Arctic Water (blue) and Barents Sea Water (gray). For the data points, each color palette represents a sampling month. The tone variations within each color palette represent the different water layers with darker tones representing deep layers and lighter tones surface layers. oceanographic Polar Front, but shared water mass characteristics in deeper layers.

Success Rate of DNA Barcoding
We had good quality amplification of the DNA extracts, even without the cleaning steps found in most extraction methods, perhaps due to the small size of the samples. A total of 1,672 individuals was barcoded, with 1,192 (71%) resulting in successful DNA-based identification (Supplementary Data Sheet 1). The success rate of DNA barcoding was very high for the majority of taxa: Echinodermata (94%), Bivalvia (93%), Bryozoa (95%), and Nemertea (93%). Annelida had 87% success and the unsuccessful 13% were caused by insufficient reads (6%), a sequence with no match in the DNA databases (3%) or by a DNA-based identification which did not match the visual identification (4%). Gastropoda had the lowest success rate (59%), most often due to the insufficient number of DNA reads in the sequencing process (35%). Approximately 86% of MOTUs (molecular operational taxonomic units) had an "exact" match (>98%) in BOLD or GenBank (Supplementary Table 1). Of those "exact" matches, 87% of unique sequences were identified to species, 8% to genus, 3% to family, and 2% to order level.

Meroplankton -Overall Diversity
Swarm 2.0 clustered the barcoded individuals into 86 different meroplankton MOTUs, assigned to 72 taxa, belonging to eight phyla (Supplementary Table 1). Taxon richness ranged from 5 to 27 taxa per sampling event (integrated over the whole water column). The species accumulation curve (Figure 3) did not reach an asymptote, indicating that we did not capture the full diversity of the system. Mollusca, Echinodermata, and Annelida dominated the meroplankton community in terms of taxon richness, but we also found some Nemertea, Bryozoa,  Table 1 and Figure 4).

Meroplankton Seasonality -Abundance
Meroplankton abundance fluctuated seasonally but patterns varied across taxa. Meroplankton was present in the water column year-round but, contrary to our initial hypothesis, total abundance did not peak around the spring bloom (Figure 4). In fact, for the location north of the Polar Front, abundance was at its lowest in June, with values around 150 individuals/m 2 . Abundance increased in August and peaked in November when it was estimated at 14,671 individuals/m 2 . Abundance was low in January (1,369 individuals/m 2 ) and reached a secondary peak in April reaching about half of November values (5,736 individuals/m 2 ). At the location south of the Polar Front, overall abundance was at its lowest in January (3,149 individuals/m 2 ) and gradually increased through the year, reaching its highest value in August (391,137 individuals/m 2 , no samples were taken in November). Total meroplankton abundance also differed on either side of the Polar Front with overall much greater abundance to the south (3,149-391,137 individuals/m 2 ) compared to the north (150-14,671 individuals/m 2 ), though this difference was not statistically significant (p = 0.206), perhaps due to the low power of the test.
These overall seasonal patterns were in large part driven by molluscs (mostly bivalves) which dominated abundance in most months, except for April when echinoderms made up most of the larval abundance, both north and south of the Polar Front ( Figure 4A). Echinoderms and annelids were also abundant and followed a similar seasonal pattern to bivalves, except that annelids peaked in June at the location south of the Polar Front and echinoderms peaked in April at the location north of the Polar Front. Gastropod larvae peaked in the winter months, in November to the north and in January to the south of the Polar Front (Figure 5). Nemerteans were present all months except April north of the Polar Front but only in August to the south. Sipunculids were only found north of the Polar Front, in November and April. Bryozoans, cnidarians, and arthropods were mostly present in January but only in low numbers.

Meroplankton Seasonality -Taxon Richness
Taxon richness was significantly greater to the north (p = 0.002). Out of the 72 taxa identified, 27 taxa were shared between the two locations whereas 20 were only found to the south and 25 only found to the north of the Polar Front (Figure 6 and Supplementary Table 1). Taxon richness also changed significantly through time (p = 0.013). North of the Polar Front, species richness was highest in November (27 taxa), thereafter decreasing gradually until reaching its minimum in June of the following year (five taxa, Figure 4B). The pattern appeared more erratic to the south, with overall high taxon richness throughout the year (>18 taxa each month) except in April when richness was lower (9 taxa).
Taxon richness was dominated by annelids with 28 taxa represented, all of which polychaete worms (Supplementary Table 1). During most months, only 2-4 polychaete taxa were found at each location, but a large rise in diversity occurred in November in the north and in June in the south when 11 and 10 taxa were identified, respectively. As with abundance, molluscs made up a large proportion of species richness at both locations at most sampling events. While bivalves made up the bulk of the mollusc abundance, however, gastropods were much more diverse with 18 different taxa identified (compared to seven bivalve taxa, Supplementary Table 1), almost all of which were nudibranchs. The bivalve Hiatella sp. K made up a large proportion of Mollusca abundance in all months, except in April when it was not observed (Figure 5A). Only one sea urchin larva (Strongylocentrotus pallidus), one sea star (Ctenodiscus crispatus), and one sea cucumber (Molpadia borealis) were found so that brittle star larvae made up the bulk of echinoderm abundance and diversity (11 species total, Figure 5C). The brittle star Ophiocten gracilis was extremely abundant and made up >80% of echinoderm abundance during all months when ophiuroids were found south of the Polar Front. It was also present to the north, but its relative contribution to the echinoderm abundance was lower there.

Larval Duration
Most taxa were present at a single time point, but a few species were found during most of the year (Figures 5, 6). The bivalve Hiatella sp. K showed no distinct trend in size through time, suggesting continuous reproduction (Figure 7). The sizefrequency histogram for the brittle star O. gracilis, on the other hand, showed a gradual increase in size through time from April to August, suggestive of a long-lived planktonic phase originating from a discrete spawning event. The bivalve Mya truncata, and the ophiuroids Ophiopholis aculeata and Ophiocten sericeum were also present in the water column throughout most of the year (Figures 5, 6), but we had insufficient numbers to generate informative size histograms.

Meroplankton Community Patterns
The final CCA model was significant (p < 0.01) and included only sampling month as a constraining factor, suggesting that seasonal differences overpowered any variation between locations or water masses (Supplementary Table 2). No strong differences between North/South locations were observed on the ordination, though visually, there seemed to be a North/South difference within months in January and June (Figure 8). The CCA explained 34.2% of total inertia in the data, with the first two axes accounting for 22%. The samples collected in November, January, and April each formed their own distinct groups on the ordination, while the communities in June and August were very similar in species composition (Figure 8). In November, the community was characterized by the polychaetes Sphaerodorum flavum, Nothria conchylega, and Phyllodoce groenlandica as well as Nudibranchia indet, the ophiuroid Ophiacantha bidentata and the bivalve Serripes groenlandicus 94.1%. The community in January was almost entirely characterized by nudibranchs. Nudibranchs were still important in April in addition to the seastar C. crispatus and the polychaetes Spio sp. and Lumbrineris sp. During June/August the bivalves Macoma calcarea, Ciliatocardium ciliatum, Mya truncata, and Hiatella sp. K as well as the ophiuroids Ophiura robusta, Ophiopholis aculeata, and O. gracilis and the polychaete Galathowenia oculata characterized the meroplankton community.

Success of DNA Barcoding
We implemented a simple, relatively cheap and highly effective method (Ershova et al., 2019) to efficiently sequence hundreds of DNA-barcodes of meroplankton individuals from eight phyla. The HotShot extraction method (Truett et al., 2000) has proven to be extremely effective in extracting DNA from meroplankton, yielding higher success rates (71% overall but >93% for Echinodermata, Bivalvia, Nemertea, and Bryozoa) than earlier studies using different extraction methods [35% in Heimeier et al. (2010); 49% in Brandner et al. (2017); 68% in Walczynska et al. (2019)]. The benefits of the HotShot extraction method include that no membrane binding of DNA and no transfer of extract are needed, and it can be scaled to high-throughput (Meissner et al., 2013). Our lower success rate for DNA barcoding of gastropods compared to the other taxa may be due to incomplete crushing of the shell or the fact that most gastropod larvae were very small (∼100 µm) and may have provided insufficient material for DNA extraction. In addition, the majority of these were suspected to be young Limacina helicina (a pelagic gastropod, not meroplankton) based on the morphological resemblance to successfully-barcoded individuals of the same species, and these can be especially difficult to DNA-barcode (Elizaveta Ershova, IMR, personal communication). Therefore, the success rate FIGURE 5 | Composition of the four most abundant meroplanktonic higher taxa: (A) Bivalvia, (B) Gastropoda, (C) Echinodermata, and (D) Annelida expressed as proportional abundance (colored bars, left y-axis) as well as total abundance (black line, right y-axis) for the locations north and south of the Polar Front. For Gastropoda and Annelida, taxa that did not make up at least 10% of the group's abundance for at least one sampling point were combined in a "Other" category for ease of reading. Some taxa are named with letters or numbers instead of a species name (e.g., Hiatella sp. K and N. conchylega CMC02). These represent currently un-named but genetically-distinct lineages.
of DNA barcoding of strictly-meroplanktonic gastropods was probably much higher than reported here.
In addition to effective DNA extraction, amplification, and sequencing, successful identification of larvae through DNA barcoding requires a curated barcode database such as BOLD against which to compare the larval sequences (Hebert et al., 2003). Important efforts in recent years (e.g., Layton et al., 2014Layton et al., , 2016 have added many additional species to the databases, but marine invertebrates in general remain underrepresented (Weigand et al., 2019). In our study, 14% of all FIGURE 6 | Heatmap of meroplankton species abundance integrated over the entire water column (individuals/m 2 ) north and south of the Polar Front. For readability, only taxa that contributed at least 1% of total abundance of at least one sampling event are included (for all taxa, see Supplementary Table 1). Taxa are clustered according to their seasonal and geographical distribution with colors indicating phylum. Note that no samples were taken south of the Polar Front in November 2017.  The blue month labels are centroids for the corresponding months. The first axis explains 12.3% of the variance while the second explains 9.6%. Frontiers in Marine Science | www.frontiersin.org sequences did not find an "exact" match (>98% sequence match) in BOLD or GenBank. Of the sequences that did get a match >98%, about a tenth was not identified to species level in the barcode database itself. Therefore, the addition of more boreal-Arctic and Arctic species to the DNA barcode databases as well as taxonomic expertise for proper identification of these organisms would further improve our ability to identify meroplankton. Our very high success rate of meroplankton identification clearly demonstrates, however, that DNA barcoding of the COI gene marker is already a highly effective way of identifying meroplankton and will continue to improve as barcode databases become more complete.

Meroplankton Biodiversity
DNA barcoding enabled identification of meroplankton at a much higher taxonomic resolution than was achieved through morphological identification. We therefore captured a greater taxon richness than reported in most other seasonal high-latitude studies using morphological identification (Schlüter and Rachor, 2001;Fetzer and Arntz, 2008;Kuklinski et al., 2013;Silberberger et al., 2016;Michelsen et al., 2017). Bowden et al. (2009) found an even higher diversity of meroplankton taxa in the Antarctic, but they sampled more frequently throughout the year and had more stations compared to the present study, possibly capturing a larger proportion of the overall meroplankton community.
In our study, molluscs, echinoderms, and annelids dominated the meroplankton community in terms of taxon richness. While annelids and echinoderms are often amongst the most diverse groups in high-latitude meroplankton assemblages (Schlüter and Rachor, 2001;Michelsen et al., 2017), to our knowledge our study is the first to identify such a diverse larval nudibranch fauna. Echinoderms, molluscs, and annelids also dominated the meroplankton community in abundance, similarly to other studies in high-latitude regions (Schlüter and Rachor, 2001;Sewell and Jury, 2011). These phyla are also abundant and taxonrich on the seafloor in the Barents Sea (Carroll et al., 2008;Cochrane et al., 2009;Jørgensen et al., 2015). Other taxa that are abundant on the seafloor in our study region (Zakharov et al., 2020), however, were notably rare or absent in our samples. First, species without a pelagic early life stage or those which remain in the benthic boundary layer just above the seafloor would not have been captured by our sampling method. A majority of polar prosobranchs, for instance, rely on direct development (Clarke, 1992) explaining why so few non-nudibranch gastropods were caught here, despite being present on the seafloor (Jørgensen et al., 2015). Additionally, our seasonal sampling may have missed species with a short-lived planktonic stage of less than a few weeks (Shanks, 2009) such as the tunicate Styela rustica (Khalaman et al., 2008). Large and mobile taxa such as crab larvae may have been able to avoid the slow-moving fine-mesh nets used in this study (Porter et al., 2008) or been too rare to be captured in a single net tow. Finally, our offshore location may explain the relative paucity of coastally dominant cirripeds [but see Basedow et al. (2010)].
While the majority of taxa found in this study are known inhabitants of the Barents Sea (Sirenko et al., 2001), a few species have, to our knowledge, never been recorded there.
Notably, the nudibranch Bohuslania matsmichaeli (99.68% DNA match) is currently only known from one fjord in southern Norway (Korshunova et al., 2018). The presence of its larva in the Barents Sea suggests either that its propagules can disperse over very large distances or that the geographical range of the adults extends further north than currently known, although a third possibility of transport by ship ballast water cannot be excluded. Note, however, that our observation of this species here is only based on one barcoded larva. The brittle star O. gracilis (100% DNA match), the larvae of which were found in this study in very large numbers, has also not been recorded in the northern Barents Sea, but is commonly found on the shelf and continental slope along the coast of Norway [Sirenko et al., 2001; Global Biodiversity Information Facility (GBIF), 2019b]. While its primary distribution is along bathyal depths (Gage and Tyler, 1981), several adults have been collected from shallower locations [Global Biodiversity Information Facility (GBIF), 2019b] so the species' ability to settle and survive on the Barents Sea shelf in the future seems unlikely but possible. It is also possible that adults of O. gracilis are currently found on the northern Barents Sea shelf but have been confused for the closely-related O. sericeum. Finally, we barcoded several individuals of the nemertean Cephalothrix iwatai (99% DNA match), a deep-sea species previously only recorded in the Sea of Japan [Chernyshev, 2013;Global Biodiversity Information Facility (GBIF), 2019a]. It is, however, highly unlikely that the larvae have drifted into the Barents Sea from the Sea of Japan. Therefore, one can presume that C. iwatai is more widely distributed than previously thought or that our specimens belong to a closely-related species which has not yet been barcoded. So far, conditions on the Barents Sea shelf have likely prevented some of these and other more southern species from settling onto the seafloor and growing to a reproductive stage, at least not in sufficient numbers to grow a detectable population. However, climate change and Atlantification of the Barents Sea (Polyakov et al., 2017(Polyakov et al., , 2020 may increase the flow of Atlantic Water onto the Barents Sea shelf , potentially transporting more propagules north as well as making conditions more favorable for some nonnative species to settle and grow there as has been suggested for other species (e.g., Berge et al., 2005;Geoffroy et al., 2018;Renaud et al., 2019).

Seasonality
Larval abundance and composition varied sharply throughout the year. Contrary to our initial hypothesis, larval abundance peaked in August (south) and November (north) rather than in the spring when food would have presumably been most abundant for the larvae. There was, however, a secondary peak in larval abundance in April north of the Polar Front. In 2018 in our study area, the phytoplankton bloom peaked in early May, first north of the Polar Front, then, about 1 week later, to the south (Marie Porter, SAMS, personal communication). Thus, our April sampling captured the beginning and our June sampling captured the end of the bloom. A mooring study at 81 • N on the northern coast of Svalbard similarly found the bloom to begin in early May . At first glance, this asynchrony of the bulk of the meroplanktonic community with the peak in primary production appears to contrast with studies in Svalbard fjords where overall meroplankton abundance peaks in spring/summer around the phytoplankton bloom (Kuklinski et al., 2013;Stübner et al., 2016), as well as in fjords in mainland Norway (Michelsen et al., 2017), in Greenland (Arendt et al., 2013), in the Antarctic and sub-Antarctic (Bowden et al., 2009;Presta et al., 2019), and off of the coast of the United Kingdom (Highfield et al., 2010). However, seasonality differs across taxa and a peak in abundance around the phytoplankton bloom is often driven by a few highly abundant taxa such as cirriped larvae (Węsławski et al., 1988;Highfield et al., 2010), which were not captured in large numbers in our study. Indeed, larvae of molluscs, echinoderms, polychaetes, and other taxa are often found later in the year than cirripeds (Węsławski et al., 1988;Highfield et al., 2010), although generally not as late as the August/November peak found here. Consistent with our findings of a late summer/fall peak (August/November) in abundance, a sediment trap study in the Chukchi Sea showed a peak in meroplankton export to the seafloor in September/October (Lalande et al., 2020). Additionally, a seasonal study of meroplankton on the continental shelf off of mainland Norway showed a peak in total abundance in July to September (Silberberger et al., 2016), consistent with our peak in August and November.
A combination of factors may explain the mismatch between the meroplankton peak and the phytoplankton bloom. First, for some of the species we collected, the pelagic larva can rely on a maternally-derived yolk sac for nutrition so it does not need to eat while in the plankton (e.g., lecithotrophic larvae of Capitella capitata, Méndez et al., 2000). For these larvae, availability of food in the plankton is irrelevant so there would be no need to time larval release with the phytoplankton bloom, at least not for the sake of the larval stage survival itself. Most taxa that we collected, however, are planktotrophic and would presumably require food to survive while in the plankton (Arctic Traits Database 4 , Degen and Faulwetter, 2019). These feeding larvae may rely on different food sources (Boidron-Métairon, 1995) than the dominant diatom taxa such as Thalassiosira spp. and Chaetoceros spp. present during the spring bloom in the study region (Hegseth, 1998;Wassmann et al., 1999). A wide range of other planktonic organisms including dinoflagellates, cercozoans, and fungi feed copepods outside of bloom periods in the Arctic (Cleary et al., 2017) and may feed meroplankton as well. Bivalve larvae, for instance, can feed on ciliates or dinoflagellates (Lindeque et al., 2015), which are present yearround in an Arctic fjord . Metabolic rate decreases with decreasing temperature in polar invertebrates (e.g., Heilmayer and Brey, 2003) so meroplankton in cold high-latitude waters may require less food to survive. Given that larval occurrence of most species in our study did not coincide with the spring bloom, it may be the adults instead which take advantage of the abundant food available during the bloom to fuel reproduction (e.g., Witte, 1996;Renaud et al., 2007b). For these income breeders, there would be a time lag between peak food availability and peak larval 4 https://www.univie.ac.at/arctictraits/ abundance (Kuklinski et al., 2013), just as we have observed in this study. In addition, some of the larvae may have been released around the spring bloom in a different location and taken several months to reach our study area. Finally, larvae released in the fall may benefit from the darker conditions to escape visual predators, though this theory remains to be tested.
Some potential biases in our sampling need evaluating. While the exact location of sampling differed slightly across months, especially for the location north of the Polar Front, we consider artifacts of this in the data minimal since the temperature and salinity profiles indicated clear differences between the north and the south of the Polar Front across all seasons. While the November sampling occurred closest to land which may partially explain the higher abundances at that time (Bluhm et al., 2018), August sampling, which also had very high meroplankton abundances, occurred furthest away from land. It is therefore likely that the August/November meroplankton peak remains valid. A secondary peak in abundance also occurred in April to the north of the Polar Front. This peak, largely driven by the ophiuroids O. gracilis and O. sericeum as well as the polychaete Spio sp., may represent a real surge in abundance just prior to the phytoplankton bloom (in May, Dalpadado et al., 2020;Dong et al., 2020). However, it could also be an artifact of the closer proximity to the Polar Front compared to other months, given that plankton often accumulate in patches around fronts (Trudnowska et al., 2016) and adult benthic invertebrates on the seafloor also occur in higher densities near the Polar Front (Carroll et al., 2008). However, O. gracilis does not occur on the northern Barents Sea shelf at all, so higher adult densities are unlikely to explain the April peak north of the Polar Front. Given that the zooplankton net mesh size was larger in April (180 µm) and, at the northern location, in June (180 and 64/300 µm, see Table 1 note for details), smaller larvae (e.g., earliest stages of Hiatella arctica measure around 120 µm, Flyachinskaya and Lesin, 2006;Brandner et al., 2017) could have been undersampled and abundances under-estimated compared to other months. Nonetheless, it is unlikely that an abundant species would have been missed entirely as some individuals would have presumably been caught in the net through aggregations of other plankton such as phytoplankton clusters. Finally, zooplankton distribution can be spatially patchy, with plankton concentrations sometimes varying by more than one order of magnitude over a distance of a few kilometers (e.g., Trudnowska et al., 2016) and meroplankton is no exception (Kersten et al., 2019). Taking multiple replicates of each net sample would help obtain a more robust estimate of species abundance and composition at a particular location. In this study, as is often the case during oceanographic surveys on shared research vessels, sampling time was limited, allowing for collection of only one replicate per sample. We recommend that future studies put a strong emphasis on replication of net samples.
The meroplankton community varied greatly throughout the year with most species present at one or two consecutive sampling events, suggesting a larval period of at most a few months for most species (assuming at least some local retention). However, some taxa (the ophiuroids O. gracilis, O. sericeum, and O. aculeata, the bivalves Hiatella sp. K and M. truncata, and the nudibranchs Dendronotus patricki and Dendronotus kalikal) were present through most of the year. Our data suggest that the extended presence in the water column can be explained in some taxa by a long-lived pelagic phase [e.g., O. gracilis, also observed in Gage and Tyler (1981)]. On the continental slope off of the United Kingdom, O. gracilis spawns in late winter and early spring and settlement of the post-larvae to the seafloor occurs in summer (Tyler and Gage, 1982). Similarly, we observed a pulse of small larvae in April and on average larger larvae as well as several post-larvae in August, suggestive of a similarly long larval phase in the Barents Sea compared to the coast of the United Kingdom, even though we would have expected a longer larval phase in our region as colder waters tend to extend larval period (O'Connor et al., 2007). Contrary to O. gracilis, the average size of Hiatella sp. K showed no obvious pattern over time, suggesting continuous reproduction (or at least multiple reproductive events) throughout the year, which is consistent with prior observations (Garcia et al., 2003;Brandner et al., 2017). In addition, several species identified in our study, such as the polychaetes Bylgides groenlandicus and N. conchylega as well as several nudibranchs, were identified in the meroplankton but, to our knowledge, a pelagic life stage has never been explicitly recorded in these species (Supplementary Image 1). Data from this study therefore contributes valuable, although coarse, species-level information on reproductive strategy, timing, and duration.

Role of the Polar Front
While seasonality drove the bulk of the differences in community composition among samples, there were some differences between the sampling locations north and south of the Polar Front as well. Specifically, about two thirds of the taxa were found only on one side of the front or the other and the taxon richness was significantly greater to the north, partially supporting our hypothesis of different meroplankton communities across the front. Note, however, that the difference in community structure was not statistically significant as per the multivariate analysis. The distinctiveness of zooplankton communities across fronts has been demonstrated in various systems worldwide (Perry et al., 1993;Coyle and Pinchuk, 2005;Powell et al., 2015), including in the Barents Sea where the Polar Front has been interpreted to act as a barrier to larval dispersal (e.g., Schlüter and Rachor, 2001). The benthic macrofaunal communities differ substantially across the Barents Sea Polar Front (Carroll et al., 2008), potentially due to different larval supplies or, vice versa, supplying larvae of different species on either side of the front. Indeed, in our study, early life stages of polychaetes (a predominantly macrofaunal group) differed greatly in species composition across the Polar Front. Of the 28 polychaete taxa identified in the meroplankton, only seven were shared between north and south, but with all phyla combined, about a third of the taxa were shared. This observation is congruent with the benthic megafaunal communities (including most echinoderms and some molluscs) on the seafloor below, which differ somewhat between the two locations but are overall more similar to each other compared to other regions of the Barents Sea (Jørgensen et al., 2015;Zakharov et al., 2020). Importantly, while our two locations were positioned on either side of the Polar Front, with surface waters differing greatly in temperature and salinity, deeper layers were similar between the two locations. The Barents Sea water mass, which spanned the deeper layers of both the north and south locations, potentially acted as a corridor for transport of larvae between our two locations. Other oceanographic features (particularly eddies, Porter et al., 2020) may also facilitate exchange of water, and therefore of meroplankton, across the oceanographic Polar Front. Eddies have been shown to transport larvae in other regions of the world (e.g., Madagascar, Ockhuis, 2016). The branch of Atlantic Water that wraps around the western coast of Svalbard and enters the Barents Sea from the north (Lind et al., 2012) is of relatively low volume and has limited penetration southward into the Barents Sea, so is unlikely to constitute a significant supply of meroplankton to our sampling site north of the Polar Front. Therefore, our study suggests that the Polar Front may not act as an impervious barrier to larval dispersal, yet maintains some level of community differentiation.
Despite some differences in faunal characteristics across the Polar Front, there was little impact of water mass on the meroplankton community composition. The same taxa were often found at all depths, in water layers of different properties, within a given station. While another meroplankton study in the Barents Sea showed strong associations of meroplankton communities with water masses (Schlüter and Rachor, 2001), another study in the Ross Sea found no significant differences across water masses once other factors, such as depth and region, were accounted for Gallego et al. (2014). A possible explanation could be that some taxa may actively migrate within the water column as a means of regulating their range of dispersal (Largier, 2003). Alternatively, we may have captured either the upward trajectory of newly-released larvae or the downward journey of settling propagules. It therefore appears that early life stages of some of these taxa may be able to withstand relatively large variations in temperature and salinity, at least for short periods of time (Costlow and Bookout, 1969).

CONCLUSION
Our new knowledge gained from the combination of high taxonomic and seasonal resolution opposes the previous notion of tight coupling of the spring bloom to larval abundance, and instead shows taxon-specific seasonal patterns with most species occurring in late summer and fall. This study also highlights the value of high-throughput DNA barcoding as a means of identifying meroplankton to high taxonomic resolution. In addition to providing species-specific information on reproductive timing, etc., this method enabled the identification of species that are non-native to the study area. Boreal species which already have a means of reaching the Arctic inflow shelves during their larval phase may be among the first to expand their range into the Arctic as the climate continues to warm. Regular monitoring of the meroplanktonic community on the inflow shelves using high-throughput DNA barcoding may play a key role in early detection of these newcomers.

DATA AVAILABILITY STATEMENT
The raw sequencing datasets presented in this study are publicly available in the Sequence Read Archive (SRA) repository of NCBI Bioproject name: PRJNA725248; https://www.ncbi.nlm.nih.gov/ bioproject/PRJNA725248. The rest of the data, as well as photos of meroplankton individuals and the bioinformatics and data analysis R codes can be found at Dryad Digital Repository https: //doi.org/10.5061/dryad.n8pk0p2vf.

AUTHOR CONTRIBUTIONS
RD, EE, PR, and BB conceived the study. RD, EE, OW, and KP designed the molecular protocol. FC and RD collected and interpreted the physical data. RD and EE collected the zooplankton samples. RD, EE, and OW performed the laboratory work. OW ran the bioinformatics and RD performed data analysis. RD wrote the manuscript with contributions from EE, OW, KP, PR, FC, and BB. All authors have read and approved the manuscript.