dataesr/scientific-paragraphs-categorization huggingface dataset:
This is especially true for groups like bivalves  (Garland and Zimmer 2002) , or in taxa comprising cryptic species which are numerous among marine invertebrates  (Appeltans et al. 2012 ). To overcome these issues, several single-species DNA-based tools have been used for the identification of NIS larvae (e. g. ,  Darling and Tepolt 2008; Harvey et al. 2009; Le Goff-Vitry et al. 2007; Sánchez et al. 2015). However, these approaches allow the identification (and sometimes quantification) of only one or a few target species. The DNA metabarcoding approach consists of the high-throughput sequencing (HTS) of selected barcodes obtained from environmental or bulk-DNA, and their taxonomic assignment based on a reference sequence database. It can identify simultaneously a large number of species from many specimens, thus being a promising method to study NIS larvae in the plankton  (Comtet et al.
Data processing After sequencing, reads were processed using the OBITools v1. 2. 11 pipeline  (Boyer et al. 2016). Briefly, paired reads were assembled and then grouped by replicate sample (i. e. , demultiplexed) according to their unique tag pair. The primers and tags were then removed and all sequences outside of a specified size range were discarded (350-450 bp for 18S and 300-320 bp for COI). Finally, singletons (sequences present only once and in only one of the PCR replicates) were discarded. PCR and sequencing errors were detected using obiclean with a value of 2 for the parameter -d (maximum number of differences allowed for two reads to be considered deriving from one another) and a value of 0. 025 for the parameter -r (threshold ratio between counts of two reads under which the less abundant is classified as an error deriving from the more abundant).
Taxonomic assignment Given the occurrence of numerous errors and missing data in public databases  (Briski et al. 2016; Harris 2003) , a custom-designed database was used for taxonomic assignment. This database was specific to the study area and curated with caution. A list of 670 species across the four target taxonomic groups (ascidians, bryozoans, bivalves and gastropods) was considered for inclusion in our custom-designed database. These species were native, non-indigenous or cryptogenic. The native species list came from "Fauna inventories of the Station Biologique de Roscoff" (available at http://www. sb-roscoff. fr/fr/observation/biodiversite/especes/inventaires/ inventaires-de-la-faune-et-de-la-flore-marines) and was completed with new data from local surveys and inventories. The non-indigenous and cryptogenic species list was composed of species that had already been reported in the study bay or were known to be present in the surrounding areas, based on both reports made by French authorities for the European Marine Strategy Framework Directive reporting, and local surveys and inventories (e.
g. ,  Bishop et al. 2015b). All available 18S and COI reference sequences associated with these taxa were retrieved from public databases (GenBank:  Benson et al. 2013; BOLD: Ratnasingham and Hebert 2007; SILVA: Quast et al. 2013; PR2: Guillou et al. 2013 ) and manually checked. Only good quality sequences, with no more than one ambiguous base, were included. In order to increase taxonomic coverage, incomplete sequences were also added, where a length of at least 75% of the complete sequence was present. Reference sequences produced locally were also included, in Table  1. Number of taxa (with the number of NIS among them indicated in bold and in parentheses), number of unique variants and number of reads identified at different taxonomic levels (species, genus, family or class) for each marker (18S and COI), within the four taxonomic groups of interest (Gymnolaemata, Gastropoda, Bivalvia, Ascidiacea) and using a BLAST ® approach against a custom-designed database.
The number of taxa (at the species, genus and family levels) with at least one reference sequence in the custom-designed database is given in the "Taxa with available references column", with the number of reference for NIS in bold and in parentheses. particular for taxa with missing data in the above cited public databases (see Data availability section). For a better discrimination of known or expected NIS, reference sequences for species belonging to the same genus were also included. Finally, all species names were checked with the World Register of Marine Species (WoRMS; WoRMS Editorial Board 2019) to ensure that all names used were valid. In the case of the Pacific oyster, both Magallana gigas and Crassostrea gigas are accepted names, so we chose to use Crassostrea gigas (Thunberg, 1793) as advised by  Bayne et al.
(2017). In the end, our custom-designed reference database was composed of 408 and 3131 sequences covering 314 and 410 species, for 18S and COI, respectively (Table  1 ). Our dataset included 42 different NIS, 35 with reference sequences for both markers; the remaining seven lacked a reference for 18S (Table  1 , Table  S1 ). Unique variants resulting from the OBITools pipeline were then compared to sequences from our custom-designed reference database using BLAST®  (Altschul et al. 1990 ). Only unique variants found in at least two PCR replicates of the same sample were considered for taxonomic assignment. Only alignments covering at least 99% of the reference sequence were considered. Variants were then assigned to the species whose reference sequence had the highest identity percentage.
If two or more alignments had the same identity percentage, the variant was assigned to the lowest common taxonomic level. Identity thresholds were defined for each marker in order to consider an assignment to the species level as valid. To select an appropriate threshold, we explored the barcoding gap of our target taxa for the two markers. To that end, each reference sequence was aligned with all the references in the database using the EMBOSS needle global alignment tool (version: 6. 5. 7. 0; Rice et al. 2000). For each alignment, the identity percentage was calculated using R 3. 4. 4 (R Core Team 2018) and the values of intraspecific and interspecific identity were plotted as a density curve using the R package ggplot2  (Wickham 2016). For COI a clear gap between intra-and interspecific identity was observed at 92% (Figure  S2B ), a threshold thus selected for species assignment.
For 18S, the distribution of the intraspecific identity is aggregated towards 100% with a clear gap visible at 98% (Figure  S2A ). However, to avoid false positive, we chose to apply a more conservative threshold of 99% for this marker. Note that the number of NIS detected using 18S was the same with both thresholds (not shown). To assess the value of using a custom-designed database for NIS detection, we also performed a taxonomic assignment with two other methods frequently used in the literature and based on public databases, either with a BLAST® approach or by using the ecotag command implemented in the OBITools pipeline. For the first method, unique variants were aligned against the GenBank nucleotide database (nt accessed 11. 01. 19) using BLAST® (v. 2. 7. 1+). Only alignments with at least 99% (18S) or 92% (COI) identity and 99% query cover were considered.
Each unique variant was then assigned to the taxon name of the reference sequence with the highest identity percentage. If several reference sequences had identical identity percentages, the unique variant was assigned to the lowest common taxonomic level. For the second approach, we used the ecotag command of the OBITools package  (Boyer et al. 2016 ) with a reference database created using the ecoPCR tool available in the same package by retrieving all sequences from GenBank, BOLD, Silva or PR2, and by adding all reference sequences produced locally. In the ecotag approach, unique variants are aligned with references using a global alignment algorithm, thus requiring only full-length reference sequences to be included in the database. Then, unique variants are assigned to the taxonomic level corresponding to the lowest common ancestor between all reference sequences that are closer to each other than to the selected unique variant  (Boyer et al.
Cross-validation with traditional methods In order to compare the results from the metabarcoding approach to more traditional ones, two additional datasets were used. They consisted in the parallel sampling, identification and quantification of larvae associated to the bivalve Ruditapes philippinarum (Adams and Reeve, 1850) and the gastropod Crepidula fornicata (Linnaeus, 1758) during the year 2012. For R. philippinarum larvae, a traditional barcoding approach was used. Briefly, an additional monthly plankton sample was collected at the study site, from March to November 2012, using the same sampling protocol. In the laboratory, each sample was adjusted to a volume of 150 mL with 96% ethanol. Bivalve larvae were then sorted in three 5 mL subsamples (10% of the sample) under a dissecting microscope aiming to randomly sort at least 100 bivalve larvae, when possible, depending on their overall abundance.
PCR products were then sequenced in both directions by Sanger sequencing. Overall we obtained an amplification-sequencing success of 64% of all sorted larvae. Obtained sequences were then compared using BLAST®  (Altschul et al. 1990 ) to reference sequences from GenBank and sequences produced locally for local species. Only alignments with 100% cover were considered. A larval sequence was assigned to a species if it differed by less than 2 base pairs (bp) (99. 8% identity) from reference sequence(s) of a single species. In all other cases (difference of 2 bp or more and/or several possible species), then the larval sequence was assigned to a family. This 2-bp threshold was based on  Blaxter et al. (2005)  and empirical observations. In this study, we only focused on larvae assigned to R. philippinarum, and results were expressed as their relative abundance compared to all bivalve larvae.
Out of these unique variants, 556 and 2,144 were assigned to an accepted species, genus or family within the four taxonomic groups of interest, when compared to the reference sequences in our custom-designed database, for 18S and COI, respectively (Table  1 ). These accounted for a total of 317,070 and 395,533 processed reads meaning that ca. 13% and 27% of the 18S and COI datasets, respectively, were assigned to one of the target taxa (Figure  1A ). With our method, all COI unique variants were assigned to the species level whereas some 18S unique variants were only assigned to higher levels (genus, family or class; Table  1 ). Overall, within the four taxonomic groups of interest, 86 and 79 taxa were identified at the species level, with 18S and COI respectively (Table  1 ). Most identified species belonged to Bivalvia (18S) and Gastropoda (COI) (Figure  1A, B ), which are the two classes with the highest number of species with a reference sequence in our database (Table  1 ).
Using alternative methods for taxonomic assignment (i. e. , BLAST® against the GenBank nt database or assignment with the OBITools suite) did not change the proportion of reads assigned to the four target classes for COI (ca. 25% of reads assigned) (Figure  1A ). For 18S, alternative methods allowed an increase in the proportion of reads assigned to a taxon (up to 40% with ecotag) (Figure  1A ). However, the number of species identified with 18S was similar across the methods, although the distribution Table  2. Non-indigenous species (NIS) in the four target classes detected by 18S, COI or both in at least one plankton sample collected in the bay of Morlaix. For each marker, the number of unique variants is given, with the number of reads in parentheses. For each NIS, the type of dispersal mode (Dispersal) is indicated with short and long disperser describing species with a life cycle including a pelagic larval stage lasting less or more than 2 days, respectively.
'Reported' indicates if the species has previously been reported in the study bay. For each marker, the total number (in bold) of reference sequences, retrieved from public databases or produced locally (in parenthesis), available in the custom-designed database is specified (N ref ). Individual DNA was tested for amplification failure with the COI primers (see Figure  S3  for amplification results). In case of COI amplification failure, N/A was added to the COI detection column. of assigned species per class and/or status species varied (Figure  1B ). With COI, the number of species identified was a bit higher due to the number of identified gastropod species, the other class and/or status being similar (Figure  1B ).
Identification of non-indigenous species Of the 42 NIS of interest for which reference sequences were available, 12 were detected with at least one marker using our custom-designed database (Table  2 ). Most of them belonged to ascidians (4) and bivalves (4), whereas two were gastropods and two bryozoans. Most species exhibit a pelagic phase with either long-lived (5) or short-lived (6) larvae. The only identified species with no pelagic larval stage is Crepipatella dilatata (Lamarck, 1822). Among these 12 NIS, nine were already reported in the study area and two had no previous record. In addition, Ruditapes philippinarum was not present in the Roscoff fauna inventories but has been observed recently. The 12 NIS were identified by at least one unique variant with 100% identity (99% cover) to a reference sequence, except three with 99% identity with an 18S reference sequence, namely Crassostrea gigas, Mercenaria mercenaria (Linnaeus, 1758) and Crepipatella dilatata.
Conversely, the six species with short-lived larvae and the species with direct development were represented by a small number of reads (less than 130 per species over all samples; Table  2 ). With the two alternative assignment methods, the number of detected NIS was 25% to 50% lower than the one achieved with our custom-designed database (Figure  1B ), and no extra NIS were identified. The use of ecotag brought the lowest number of NIS (6) with five species assigned with 18S (A. humilis, C. eumyota, C. fornicata, Mya arenaria Linnaeus, 1758, and R. philippinarum) and four with COI (Bugula neritina (Linnaeus, 1758), C. fornicata, M. arenaria, and R. philippinarum). As compared to ecotag, the BLAST® approach against the GenBank nt database brought the same six NIS as well as C. gigas and Mercenaria mercenaria with 18S and B.
Comparison of the metabarcoding results with other datasets The metabarcoding results were compared with data obtained through either barcoding of individual larvae (Manila clam R. philippinarum) or morphological identification of larvae (slipper limpet C. fornicata). Whatever the method considered, R. philippinarum was observed in three samples only, namely August, September and October 2012 (Figure  4 ). The 18S reads   assigned to the Manila clam R. philippinarum accounted for between 18% and 62% of the number of reads assigned to bivalves in the samples where the species occurred (Figure  4A ). This proportion was even higher when considering COI, ranging from 67% to 95% across samples (Figure  4B ). These temporal variations matched well those displayed by the number of larvae identified using traditional barcoding, even if in this case the proportion of clam larvae was lower, ranging from 16% to 30% (Figure  4C ).
In the metabarcoding datasets C. fornicata was identified in all samples from March to October 2012 with both markers (except the 28. 08. 12 and 10. 10. 12 with 18S; Figure  5A, B ), which was congruent with the observation of larvae from this species based on morphological identification over the same year (Figure  5C ). No correlation between larval counts and the number of reads was observed, neither for 18S (r = -0. 07, p = 0. 81) nor COI (r = 0. 29, p = 0. 29).
Discussion DNA metabarcoding is predicted as a promising approach to identify NIS at all steps of the introduction process (e. g. ,  Comtet et al. 2015; Darling et al. 2017; Zaiko et al. 2018) , and this study provides evidence supporting this statement. We examined zooplankton samples collected over 22 consecutive months, targeting larvae of benthic NIS. By using two DNA markers, 18S (nuclear DNA) and COI (mitochondrial DNA), and a custom-designed database, we identified 12 NIS among the four targeted taxonomic groups (Table  2 ). Our method proved to be efficient in detecting both long and short dispersers. Moreover, three of the identified species had never been recorded in the study area and could be novel introductions, previously unnoticed ancient introductions, or false positives. Finally, seasonal variations inferred from read distributions of several species were consistent with their known reproductive periods.
Zooplankton metabarcoding: an efficient tool to detect non-indigenous benthic species Among the four classes studied, 42 NIS were targeted, present either in the study bay or in neighbouring areas, notably along the French and English Western English Channel coasts. Twelve of them (29%) were identified using our approach. However, only 17 species out of 42 were previously reported in the study bay, and one of them (Ciona robusta Hoshino and Tokioka, 1967) was observed only once in 2012 and never observed since then despite regular surveys  (Bouchemousse et al. 2017). Thus, 75% of these 16 expected NIS were actually detected. Zooplankton metabarcoding thus appears as a powerful tool for NIS detection provided that multiple markers and a custom-designed database are used. Only four of the 12 NIS were identified with both COI and 18S.
As they have all been reported in the study area, we assumed that they were actually present in our samples. Thus 18S can be a useful complement to COI, but it is critical to verify its discriminating power to avoid drawing erroneous conclusions about the presence of a NIS. Our ability to identify a species not only depends on the amplification efficiency and discrimination power of the chosen markers, but also on the quality of the reference sequences. Building a manually-checked customdesigned reference database (with good quality sequences and reliable identification) comprising hundreds of species is time-consuming, and one might wonder if the reliability of the results obtained is worth the effort. Taxonomic assignment is commonly based on public databases, like the nucleotide collection of GenBank (nt), which are recognized to include mistakes (e.
g. ,  Harris 2003) , with risks of false assignments or failure to detect target NIS. One common example of such risks is the description of a new species within a species complex. This is the case of the tunicates Ciona intestinalis (Linnaeus, 1767) and C. robusta, which were reported in our study area  (Bouchemousse et al. 2016b). Formerly known as C. intestinalis type B and type A, respectively  (Zhan et al. 2015) , they were recently reclassified as C. intestinalis and C. robusta based on molecular and morphological evidences  (Gissi et al. 2017, and references therein). No reliable assignment to the species level would have been possible since some references obtained before the taxonomic revision are still falsely attributed to C. intestinalis (e. g. , AK173399. 1). More specialized and curated databases (e.
g. , BOLD, Silva) must be preferred for species identification, if references are available for the targeted taxa. Besides, the completion of the reference database is a key condition for an accurate detection of NIS  (Briski et al. 2016; Comtet et al. 2015)  and producing new reference sequences for known or potential NIS is worth the effort. In our results, three NIS were not identified when using BLAST® with GenBank nt (Crepipatella dilatata, Watersipora subtorquata (d'Orbigny 1852), and Botrylloides diegensis). No public reference was available for our targeted markers, and they were identified based on our locally-produced references. Another species, Bugula neritina, was also not detected with the 18S marker when using the nt public database. Divergent lineages often exist within accepted species  (Pante et al.
2015) , and species complexes have been uncovered in many invasive species, such as Bugula neritina  (Fehlauer-Ale et al. 2014)  or Watersipora spp. (Mackie et al. 2012). Obtaining local references is therefore important to ensure correct assignment, particularly when using a threshold value. For example, the sequences we produced locally for B. neritina diverged by 3% from the one available in GenBank (AF499749. 1), which most likely originated from a Chinese population. This difference, higher than the threshold chosen for assignment with 18S, prevented the identification of this NIS when using only public references. Finally, the alignment method used for taxonomic assignment matters, as shown here by using two alignment tools (BLAST® versus ecotag command from OBITools). ecotag required building a composite reference database composed solely of full-length sequences, which lead to a reduction in the number of available references and the number of species with references.
Read counts as a proxy for larval abundance and NIS reproductive success Detecting reads of a specific benthic NIS in a zooplankton metabarcoding dataset most likely indicates the presence of its larvae, although an alternative origin (shed cells, mucus, see above) could not be discarded. When resulting from the presence of larvae, read occurrence would be indicative of the reproductive status of populations that have most likely established in the vicinity of the sampling area. In such cases, the temporal distribution of reads can further allow to investigate the species reproductive dynamics. In the three benthic NIS with long-lived larvae and large number of reads (i. e. , M. arenaria, R. phillipinarum, C. fornicata), such distribution suggested a seasonal pattern congruent with their known reproductive cycle.
In both datasets an abundance peak is observed in May, whereas a second peak is observed in September or in August, based on the metabarcoding and morphological datasets, respectively. In addition, when comparing larval abundances of the slipper limpet C. fornicata based on morphological identification to the number of reads assigned to this species in 2012 (Figure  5 ), no correlation was found with either marker or even between both markers. This result is congruent with what have been found in other studies (e. g. ,  Elbrecht and Leese 2015; Pochon et al. 2013; Sun et al. 2015)  and compels us to nuance our data interpretation. Many elements can explain why we did not observe a correlation between read and larval counts for C. fornicata. For instance, the comparison between read abundances and individual counts can be arguable since the size of individuals will influence the quantity of DNA available after extraction, and thus their relative abundance in metabarcoding assays  (Elbrecht and Leese 2015; Elbrecht et al.
Another possible explanation is that the larvae are exported off the bay, a hypothesis suggested to explain the low proliferation of C. fornicata at this same location  (Rigal et al. 2010). Finally, the increase in the proportion of infertile triploid Pacific oyster in French production might have contributed to the low number of larvae produced by local farming. In contrast, the native flat oyster (Ostrea edulis), which only accounts for a low percentage of the produced oysters in the bay, has been identified with a high number of reads in both datasets (Figure  3 ). The seasonal pattern displayed by metabarcoding results is in agreement with the known reproductive period of this species  (Eagling et al. 2018) , which is known to spawn at temperatures lower than those required by C. gigas  (Mann 1979).