Metabarcoding using multiplexed markers increases species detection in complex zooplankton communities

Abstract Metabarcoding combines DNA barcoding with high‐throughput sequencing, often using one genetic marker to understand complex and taxonomically diverse samples. However, species‐level identification depends heavily on the choice of marker and the selected primer pair, often with a trade‐off between successful species amplification and taxonomic resolution. We present a versatile metabarcoding protocol for biomonitoring that involves the use of two barcode markers (COI and 18S) and four primer pairs in a single high‐throughput sequencing run, via sample multiplexing. We validate the protocol using a series of 24 mock zooplanktonic communities incorporating various levels of genetic variation. With the use of a single marker and single primer pair, the highest species recovery was 77%. With all three COI fragments, we detected 62%–83% of species across the mock communities, while the use of the 18S fragment alone resulted in the detection of 73%–75% of species. The species detection level was significantly improved to 89%–93% when both markers were used. Furthermore, multiplexing did not have a negative impact on the proportion of reads assigned to each species and the total number of species detected was similar to when markers were sequenced alone. Overall, our metabarcoding approach utilizing two barcode markers and multiple primer pairs per barcode improved species detection rates over a single marker/primer pair by 14% to 35%, making it an attractive and relatively cost‐effective method for biomonitoring natural zooplankton communities. We strongly recommend combining evolutionary independent markers and, when necessary, multiple primer pairs per marker to increase species detection (i.e., reduce false negatives) in metabarcoding studies.

tion and taxonomic resolution. We present a versatile metabarcoding protocol for biomonitoring that involves the use of two barcode markers (COI and 18S) and four primer pairs in a single high-throughput sequencing run, via sample multiplexing. We validate the protocol using a series of 24 mock zooplanktonic communities incorporating various levels of genetic variation. With the use of a single marker and single primer pair, the highest species recovery was 77%. With all three COI fragments, we detected 62%-83% of species across the mock communities, while the use of the 18S fragment alone resulted in the detection of 73%-75% of species. The species detection level was significantly improved to 89%-93% when both markers were used. Furthermore, multiplexing did not have a negative impact on the proportion of reads assigned to each species and the total number of species detected was similar to when markers were sequenced alone. Overall, our metabarcoding approach utilizing two barcode markers and multiple primer pairs per barcode improved species detection rates over a single marker/primer pair by 14% to 35%, making it an attractive and relatively cost-effective method for biomonitoring natural zooplankton communities. We strongly recommend combining evolutionary independent markers and, when necessary, multiple primer pairs per marker to increase species detection (i.e., reduce false negatives) in metabarcoding studies.

K E Y W O R D S
18S, cytochrome c oxidase subunit I, metabarcoding, multigene, multiple primer pairs, zooplankton 2003; Ratnasingham & Hebert, 2007), with high-throughput sequencing (HTS) technologies to reveal species composition in "bulk" samples or environmental DNA (eDNA) samples (i.e., DNA that leaks into the environment; reviewed in Taberlet et al., 2012). Although metabarcoding is a very promising method, its efficient application is still hindered by several technical limitations which are often responsible for generating both false negatives (species being present in a sample but not detected) and false positives (species being detected but not present). This method relies on well-designed primers to amplify a homologous marker gene from a taxonomically complex sample (Creer et al., 2016). Thus, challenges often include finding a suitable DNA region to amplify across target taxa, dealing with PCR amplification errors and sequencing artifacts, developing highquality reference sequence databases, and choosing the appropriate bioinformatic steps to accommodate variable sequence divergence thresholds among species (Cristescu, 2014;Taberlet et al., 2012;Yoccoz, 2012). Choosing one or more appropriate genetic markers for metabarcoding is considered essential for accurate molecular species identification (Bucklin, Lindeque, Rodriguez-Ezpeleta, Albaina, & Lehtiniemi, 2016;Clarke, Beard, Swadling, & Deagle, 2017), as it affects both PCR amplification success and species-level resolution.
To allow efficient species identification, the genetic marker used must show high interspecific variation and low intraspecific variation. However, it is often difficult to strike a balance between high amplification success across diverse taxon groups and species-level resolution (Bohle & Gabaldón, 2012). Markers that undergo fast rates of evolution have discriminative taxonomic power for resolving closely related species but often lack conserved primer binding sites appropriate for amplifying broad taxonomic groups. Degenerate primers are often designed when conserved primer binding sites are not available. However, primer-template mismatches can generate imperfect primer match with some DNA templates (Pinol, Mir, Gomez-Polo, & Agust, 2015). This primer bias often distorts the biotic composition.
Most current metabarcoding projects use a single locus approach, and the most common markers are the cytochrome c oxidase subunit I (COI) for animals (Hebert et al., 2003;Leray et al., 2013), internal transcribed spacer (ITS) for fungi (Horton & Bruns, 2001;Schmidt et al., 2013), and plastid DNA (matK and rbcL) for land plants (Chase & Fay, 2009;Yoccoz, 2012). Alternative single markers are standardly used for particular taxa. For example, 12S is the most commonly used metabarcoding marker for fish (Miya et al., 2015;Valentini et al., 2016). Using a single organelle marker can occasionally cause erroneous species identification due to interspecific mitochondrial introgressions (Funk & Omland, 2003;Meyer & Paulay, 2005); therefore, the use of both uniparentally inherited organelle DNA and biparentally inherited DNA has been recommended (Taberlet et al., 2012). The mitochondrial COI gene has high resolution for species identification and relatively extensive reference sequence libraries (Ratnasingham & Hebert, 2007), but it is often difficult to amplify consistently across diverse taxonomic groups due to lack of conserved primer binding sites (Deagle, Jarman, Coissac, Pompanon, & Taberlet, 2014). It was suggested that using well-designed degenerated COI primers could reduce the COI primer bias (Elbrecht & Leese, 2017).
An alternative approach, tested here, is the use of multiple primers pairs per markers. In contrast to the mitochondrial COI gene, the nuclear 18S gene provides more conserved priming sites for greater amplification success across broad taxonomic groups, but often provides lower resolution for species identification (Bucklin et al., 2016;Hebert et al., 2003;Saccone, Giorgi, Gissi, Pesole, & Reyes, 1999). Another major disadvantage with using 18S as a metabarcoding marker is that it varies in length at V4 region across diverse species, causing sequence alignment uncertainties across broad taxa and consequently difficulties estimating divergence thresholds and implementing clustering approaches for species identification (Flynn, Brown, Chain, MacIsaac, & Cristescu, 2015;Hebert et al., 2003).
These marker-related problems led many researchers to propose the need to use multiple markers in metabarcoding studies (Bucklin et al., 2016;Chase & Fay, 2009;Drummond et al., 2015). The multimarker approach has the potential to reduce rates of false negatives and false positives. Despite these promises, a multigene approach has rarely been applied in metabarcoding studies, and comparisons of biodiversity estimates across the different markers are usually not reported (e.g., COI for metazoan and RuBiscCO for diatoms, Zaiko et al., 2015; species-specific primer pairs of COI and cytochrome b markers, Thomsen et al., 2012; chloroplast trnL and rbcL for surveying different terrestrial habitats, Yoccoz, 2012). In addition, many multimarker metabarcoding studies use a single primer pair per marker (Clarke et al., 2017;Drummond et al., 2015;Kermarrec et al., 2013). Using multiple primer pairs is expected to reduce amplification biases and increase the opportunities of detecting all targeted taxonomic groups. To fully understand the performance of a multigene metabarcoding approach, mock communities are ideal because the expected number of species is known a priori (Clarke, Soubrier, Weyrich, Weyrich, & Cooper, 2014;Elbrecht et al., 2016;Kermarrec et al., 2013). Nonetheless, there are few studies that have taken this approach (but see Clarke et al., 2014). As such, there is an urgent need for experiments that test species detection rates and taxonomic identification accuracy in mock communities using multimarker and multiprimer pair metabarcoding to test the validity of this method for biomonitoring.
In this study, we use a combination of mitochondrial (COI) and nuclear markers (18S) and multiple COI primer pairs in a single Illumina run for recovering species by metabarcoding mock communities of zooplankton. Species detection is assessed among markers and primer pairs to evaluate the benefit of multimarker and multiprimer pairs per marker. We also compare species detection rates and detection accuracies with a single-marker metabarcoding experiment in which 18S was used alone. We calibrate the multiplexing multigene approach using a series of mock communities containing single individuals per species (SIS), multiple individuals per species (MIS), and populations of single species (PSS). The resulting calibrated workflow performs better than a single marker or single primer pair approach and can be applied to assess zooplankton biodiversity in natural aquatic habitats.

| Primer testing
Preliminary primer amplification tests were conducted qualitatively on a total of 103 species using 13 COI primer pairs (COI-5P region) and one 18S primer pair (V4 region; see Supporting information Table S1 for the complete list of primers). We selected primer pairs known to provide amplification success across a wide range of taxa as well as good discriminatory power for species identification. The only 18S primer pair tested is known for its successful amplification across a broad range of zooplankton groups (Brown, Chain, Zhan, MacIsaac, & Cristescu, 2016;Zhan et al., 2013). Specimens used for primer testing were sampled from 16 major Canadian ports across four geographic regions (Atlantic coast, Pacific coast, Arctic, Great lakes; Chain et al., 2016;Brown et al., 2016) and were identified morphologically by taxonomists. A total of 103 species belonging to the phyla Rotifera, Crustacea, Mollusca, and the Subphylum Tunicata were selected and tested (see Supporting information Table S2). A subset of those species was used to assemble mock communities for metabarcoding validation (see Supporting information Table S2). PCR amplification was performed in a total volume of 12.5 μl: 0.2 μM of each forward and reverse primers, 1.25U taq DNA polymerase (GeneScript, VWR), 2 mM Mg 2+ , 0.2 μM dNTP, and 2 μl of genomic DNA. The PCR conditions of each primer pair were based on their sources in the literature (Supporting information Table S1). After the broad primer testing, four primer pairs were selected for metabarcoding several mock communities: one targeting the nuclear 18S V4 region (Zhan et al., 2013) and three COI primer pairs producing three different (partially overlapping) fragments within the COI-5P region ( Figure 1, Supporting information Figure S1, Table 1).

| Assemblage of mock communities
Mock communities were constructed with the aim of incorporating various levels of genetic variation (intragenomic, intraspecific, interspecific) and representing natural zooplankton communities including species from broad taxonomic groups: Mollusca, Rotifera, Tunicata, and Crustacea (Amphipoda, Anostraca, Cladocera, Cirripedia, Copepoda, Decapoda) (see Supporting information Table   S3 for species list). Morphologically identified specimens are at the species or genus level, with a few exceptions that were identified only to the family level. DNA was extracted from each specimen using Qiagen DNeasy Blood & Tissue kits and stored in ultrapure water in the freezer at −20°C as described in Brown, Chain, Crease, MacIsaac, and Cristescu (2015). The DNA was eventually combined into several different mock community assemblies. Laboratory blanks were used consistently during DNA extractions to assure no interference from contamination.   Table S3: 3a1-d3), respectively. The inclusion of single individuals in the SIS communities allowed examination of species detection with only interspecific variation. The MIS communities, which most closely resembled natural communities, allowed the examination of species detection with both intraspecific and interspecific variation.
The PSS communities allowed the examination of intraspecific F I G U R E 1 The amplified fragments used for metabarcoding. The 5′ end fragment of 325 bp refers to the FC fragment matching the COI-5P gene before the nucleotide position 400. The 3′ end fragment of 313 bp refers to the Leray fragment matching the COI-5P gene after nucleotide position 300, and the whole COI-5P gene of 658 bp refers to the Folmer fragment with forward reads R1 matching before nucleotide position 300 and the reverse reads R2 matching after nucleotide position 400. The primers are not included in the fragment lengths, and the gray lines refer to the forward and reverse reads from the paired-end 300 bp Illumina MiSeq next-generation sequencing. *The 18S fragment sizes vary between species, resulting in some forward and reverse reads that do not overlap variation and the relationship between read abundance and the number of individuals of the same species.

| Library preparation and next-generation sequencing (NGS)
DNA extractions were first quantified using PicoGreen (Quant-iT ™ Picogreen dsDNA Assay Kit, Thermo Fisher Scientific Inc.), then diluted to 5 ng/μl. The protocol "16S Metagenomic Sequencing Library Preparation" (Illumina Inc.) was used with small modifications to prepare the sequencing-ready libraries. Library preparation involved a first PCR, followed by a first cleaning with Agencourt AMPure beads (Beckman Coulter Life Sciences Inc.), a second PCR with Nextera Index kit (V3), and a second clean-up prior to next-generation sequencing (NGS). The first PCR was conducted in two replicates for each library and each of the four DNA fragments. Negative controls were included in each round of PCRs.
PCR amplification was performed in a total volume of 12.5 μl: 0.2 μM of each forward and reverse primers, 6 μl of 2xKAPA HiFi HotStart ReadyMix (KAPA Biosystems Inc., USA), and 1.5 μl of diluted genomic DNA. Due to the incompatibility of KAPA kit with primers involving inosine ("I") in the COI primer Ill_C_R (Shokralla et al., 2015), all the FC fragments were amplified using a standard PCR gradient with taq DNA polymerase (GeneScript, VWR) as in the original paper (Shokralla et al., 2015). The PCR thermocycler regimes were the same as in the original papers: 18S V4 (Zhan et al., 2013), FC (Shokralla et al., 2015), Leray (Leray et al., 2013), and Folmer (Folmer, Black, Hoeh, Lutz, & Vrijenhoek, 1994) (see Figure 1 for details). The two replicates of each PCR reaction for each fragment were pooled together after visualization on a 1% electrophoresis gel. The PCR products were quantified and pooled (equal volumes from each fragment) so that each library contained all four fragments. After this step, there was a total of 24 libraries each with four different PCR amplicons (Supporting information Table S3): six SIS libraries (simple communities); six MIS libraries (complex communities); and 12 PPS libraries (single species communities). The 24 libraries obtained were cleaned using ultrapure beads at ratio of 0.875 (28 μl beads in 32 μl solution), indexed using Nextera ® XT Index Kit (24-index, V3), and final clean-up using ultrapure beads to become sequencing-ready. All libraries were submitted to Genome Quebec for final quantification, normalization, and pooling and were sequenced using pair-end 300-bp reads on an Illumina MiSeq sequencer in one run. Note that the four single individuals per species (SIS) libraries (1a, 1b, 1c, 1d) and the four multiple individuals per species (MIS) libraries (2a, 2b, 2c, 2d) were quantified and pooled in equal molar for next-generation sequencing. The PSS libraries (3a1-d3) were quantified and pooled in different molars relative to their number of individuals of the species. It is worth nothing that the pooling of PCR amplicons prior to indexing makes this a more cost-effective approach than methods that separately index each PCR amplicon prior to pooling.
The same genomic DNA of the four SIS (libraries 1a, 1b, 1c, 1d) and the four MIS (libraries 2a, 2b, 2c, 2d) communities was TA B L E 1 The four primer pairs used in this metabarcoding study: 18S primer pair amplifying the V4 region and three COI primer pairs amplifying different fragments of the COI-5P gene.
See Supporting information Table S1 for the complete list of primers used for the preliminary primer testing step sequenced using only the 18S primers in a separate run. The library preparation and sequencing was performed in the same proportions of 5% of one flowcell using the same Illumina MiSeq pair-end 300-bp platform. This experiment was conducted to compare sequencing depth and species detection rates between a metabarcoding run with a single marker versus a multiplexed metabarcoding run with other markers/fragments (more than one marker/fragment per run for the same sample/library).

| Building a local reference database
We created a local database composed of 149 total sequences used for taxonomic assignment of reads (see Supporting information  Figure 1 for the detailed fragments positions). The 18S reference sequences contained the V4 region without trimming. The best BLAST hit against our local reference database was used to classify each sequence read with a minimum of 95% identity and an alignment length of at least 150 bp in forward and reverse reads. These relatively relaxed thresholds were used to accommodate the species with congeneric or confamiliar reference sequences. In the case of multiple best hits, if the correct species assignments could not be confirmed manually based on reads blasting against the online database on GenBank, they were excluded from further analysis.

| Bioinformatics analyses
The bioinformatic pipeline in this study consisted of demultiplexing, quality filtering, trimming raw reads, and assigning taxonomy via BLASTN (Altschul, Gish, Miller, Myers, & Lipman, 1990) against our local reference database. Taxonomic assignment was conducted at a minimum of 95% identity. We performed first a local BLAST to find unique best hits. When multiple species had equal identity, a second BLAST search in GenBank was performed to find unique best hits. If multiple species still appeared as having equal identity, the read was excluded from further analysis (Supporting information Figure S2). Each mock community was processed as a separate "library" and could be demultiplexed via their unique combination of indices. Raw reads were assigned to their corresponding libraries, generating paired forward R1 and reverse R2 files for each library (Raw read pairs in Table 2). The raw reads were then quality filtered and trimmed via "Quality Trimmer" from the FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/), with a minimum Phred quality score of 20 and a minimum length of 150 bp after trimming (see trimmed-R1 and trimmed-R2 in Table 2). After quality trimming, the R1 and R2 reads were separately used as queries in BLAST against the local database. The BLAST results of R1 and R2 were concatenated together (see paired reads after trimming in Table 2), and only the sequences with both R1 and R2 returning an accepted BLAST match to a reference sequence (>95% identity and >150 bp) were kept for further analysis (see filter-blasting step in Table 2). The BLAST results were then further filtered based on whether both R1 and R2 reads were assigned to the same species (see filter-blasting same species in Table 2 and Supporting information Figure S2).

| Primer testing
The preliminary amplification success of 14 primer pairs was tested on 103 species (Supporting information Table S2). The highest success of 76% was observed for the 18S fragment (Zhan et al., 2013), followed by the COI_Radulovici fragment (58%) (Radulovici, Bernard, & Dufresne, 2009) and then the COI_FC fragment (50%) (Shokralla et al., 2015). The overall amplification success rate of the COI_Leray fragment (38.5%) (Leray et al., 2013) was similar to the COI_Folmer fragment (37.5%) (Folmer et al., 1994). Although the three COI fragments COI_FC, COI_Leray, and COI_Folmer were designed to target a wide range of phyla, amplification success was found to be dependent on the species group. When selecting the primer for the metabarcoding study, we considered not only the overall performance of the primers under our specific conditions but also the size of the amplicons, as well as the general use of the primer pair in other barcoding-related studies.

| Read abundance comparison
A total of 20.73 million raw read pairs were sequenced from the mock communities, of which 16.72 million paired reads remained after quality filtering (Table 2)

| The performance of the 18S marker when used alone vs. with other markers
The method tested here is a multimarker approach with more than one marker sequenced in one run rather than requiring multiple runs, making it versatile and cost-effective. However, the impact of this "multiplex" approach on species detection rates and sequencing depth (number of reads per individual/species) needs to be examined. We compared results from the 18S marker in our multiplexed multimarker approach to those in a single marker approach using both the SIS communities (Supporting information Table S5) and MIS communities (Table 3). In both cases, the sequencing depth (number of reads) on average and per individual or per species was consistent between the single-marker and multimarker datasets (Table 3 and Supporting information Table S5). In the SIS communities of 56 species, discrepancies were only found when read counts were very low.
For example, three species were detected exclusively in the singlemarker dataset, while three species were detected exclusively in the multimarker datasets, but the number of reads in all six of these instances was low (≤11 reads), representing about 0.003% of the total taxonomically-assigned reads. In the MIS communities of 14 species, only two species had different detection between the single-marker and multimarker datasets: Leptodora kindtii was only detected in the single-marker experiment with 47 reads, and Corbicula fluminea was only detected in the multimarker datasets with nine reads (see Table 3). This demonstrates that the majority of species were consistently detected in both single-marker and multimarker metabar-

| Primer pair choice and species detection
The three different COI fragments selected correspond to the COI-5P region (Figure 1), encompass regions with different levels of genetic variation across the species included in the mock communities, and show variation in the amplification success of various taxonomic groups. The number of species detected among the three COI fragments was compared in both SIS communities and MIS communities (Figure 4). We found that few species (2%-3%) were uniquely recovered by a single COI fragment, and most of the species (45%-60%) were consistently detected by all three COI fragments. The three COI fragments together detected 59% of the species in SIS and 80% of the species in MIS communities ( Figure 4).
The use of three COI primer pairs improved species detection rates by 3.8%-7.5% in SIS and 3%-17% in MIS communities compared to using a single COI primer pair.

| Marker choice and species detection
Species detection rates in the SIS communities and MIS communities were also compared between the 18S V4 marker and the three COI fragments considered together ( Figure 5) (Figure 6), improving species detection rates by 11.3%-16.6% compared to using the 18S marker alone and by 13.3%-30.2% compared to using the three fragments of the COI marker (see Supporting information Table S6 for detailed species detection).

| False positives: detection of species not intentionally included in mock communities
The incidence of false positives (species detected but not intentionally included) in the three main communities were compared between the 18S V4 marker and the three COI fragments (Table 4).
In general, a low number of reads (sometime single reads) were assigned as contaminants (Supporting information Tables S7 and S8). Notes. The "n" refers to the number of individuals per species. Two instances when species detection differed between the 18S marker alone and 18S marker multiplexed are marked as an asterisk. Pearson's correlation coefficient r = 0.873, R 2 = 0.763.
The use of multiple group-specific COI primer pairs has been suggested as an efficient method for obtaining higher amplification success when studying broad taxonomic groups (Bucklin et al., 2016;Cristescu, 2014). Moreover, the use of both uniparentally inherited markers such as COI and biparentally inherited markers such as 18S has been suggested as an efficient method for increasing the accuracy of species identification (Taberlet et al., 2012). Through the use of mock communities with known taxonomic composition, we demonstrate that a multigene (COI and 18S) and multiprimer pair (three COI primer pairs) metabarcoding approach can improve species detection and provides the built-in ability to cross-validate results.

| Multiple primer pairs
The mitochondrial COI marker has been reported to be technically challenging for amplification of broad taxonomic groups due to the lack of conserved priming sites (Bucklin et al., 2016;Deagle et al., 2014). Both group-specific (Bucklin et al., 2010) and species-specific (Thomsen et al., 2012) primer pairs have been used in COI barcoding and metabarcoding. The 18S primer pair used in this study targets the V4 region of zooplankton and was successful in previous metabarcoding studies Chain et al., 2016;Zhan et al., 2013). The 13 COI primer pairs tested here showed major differences in overall amplification success depending on the group of species. Overall, amplification success of the 13 COI primer pairs followed species-specific rather than group-specific patterns in the majority of taxa tested here (Supporting information Table S2). In addition to amplification success across taxa of interest, amplicon length is also an important consideration for studies using degraded environmental DNA, which require short amplicons (Meusnier et al., 2008), and is upwardly limited by the capacity of NGS technology to obtain accurate long reads (Shaw, Weyrich, & Cooper, 2017). For example, primer pairs used here that amplified more than 600 bp (Tables S1 and S2)  Most metabarcoding studies use a single primer pair, but multiple primer pairs (species-specific or not) has been suggested and shown to improve amplification success from community samples (Bucklin et al., 2010(Bucklin et al., , 2016Clarke et al., 2014;Elbrecht & Leese, 2017;Thomsen et al., 2012). Species detection rates of the three COI fragments in our metabarcoded mock communities were expected to be higher than species amplification success during the qualitative primer testing due to massive parallel sequencing and high level of sensitivity. This was generally true but species detection varied across the three COI fragments presumably due to primer biases (see Supporting information Figure S1 for comparison).  (Clarke et al., 2014;Letendu et al., 2014;Thomsen et al., 2012). The multiple COI primer pairs covering different regions of the same marker in this study were found to improve species detection rates in both SIS (3.8%-7.5%) and MIS (3.3%-16.7%) mock communities. However, degenerate COI primer pairs have been shown to have better species detection rates than nondegenerate primers (Elbrecht & Leese, 2017) when very broad taxonomic groups are investigated. Therefore, the use of degenerate reverse primer for the Leray and Folmer fragments may farther improve the species recovery rates. The use of multiple primers pairs can be applied as an alternative approach for the markers without such fully degenerated primers available.

| Marker choice
It is generally accepted that the choice of metabarcoding marker can greatly affect species estimates (Bucklin et al., 2016;Cristescu, 2014;Tang et al., 2012). Nevertheless, only a limited number of metabarcoding studies have used a multigene approach, and the use of multiple evolutionarily independent markers has even more rarely been sequenced in a single NGS run. A few metabarcoding biodiversity studies have compared 18S and COI markers, with results varying across different taxonomic groups. Drummond et al. (2015) reported both COI and 18S markers providing good proxies to a traditional biodiversity survey dataset for soil eDNA. Tang et al. (2012) reported that COI in eDNA surveys of meiofauna estimated more species than morphospecies (species identified by morphology),

Species recovery rates
Single individuals per species (n = 53) Multiple individuals per species (n = 30) related species. Overall, the combination of 18S and COI improved species detection rates by 11%-30% compared to using a single 18S or COI marker with the tested primers.
Sequencing depth is often of major concern for fully describing community members from a complex sample. The number of libraries pooled in one sequencing run affects the number of reads per species (Letendu et al., 2014;Shaw et al., 2017). As expected, we found that the number of reads per individual or species varied significantly across markers and fragments. We consider that efficient equimolar quantifications prior to pooling including amplicons of similar length and adjusted bioinformatics pipelines could potentially also counter this variation. On a more positive note, the number of reads assigned to each species and overall species detection rates were consistent whether using a single-marker or multimarker metabarcoding approach. Therefore, the sequencing depth and species detection rates were not affected using multiple markers in one sequencing run, indicating that multiplexing several primer pairs and markers can provide a robust method to characterize samples without appreciably sacrificing read depth or species detection.
Our study compares species detection success in zooplankton metabarcoding using two evolutionarily independent markers combined with different primer pairs of the same marker. It is important to recognize that the relatively high species recovery we report might not be achieved in studies applying different bioinformatics steps such as implementing OTU clustering methods or using online reference databases which are likely to increase both false positives and false negatives. With the increasing data output from NGS technologies and the ability to pool libraries for sequencing, our results support the use of multiple genetic markers as a cost-effective approach to assessing biodiversity in a broad range of taxa within the same run. This approach also provides a built-in means to cross-validate species detection among the markers. PCR-free methods have been developed to avoid PCR bias and to enable use of more markers Zhou et al., 2013). Through this study, the use of two evolutionarily independent markers significantly improved species detection rates, and the use of three COI primer pairs improved species detection rates for particular taxa. TA B L E 4 Rates of false negatives (species not detected) and false positives (species detected but not included in the mock communities) for the four fragments in all libraries

| CON CLUS IONS
Most metabarcoding studies to date have sequenced single markers, but the choice of marker is known to greatly affect species estimates and detection accuracy. Our results suggest that a multiplexed metabarcoding approach using multiple markers and multiple primer pairs can ultimately achieve more accurate biodiversity estimates by reducing both false positives and negatives. Furthermore, the sequencing depth (number of reads per species) and species detection rates remained consistent whether multiplexing multiple fragments or using a single marker. Overall, our metabarcoding approach utilizing multiple markers and multiple primer pairs improved the species detection rates compared to using a single primer pair and/or marker. Thus, metabarcoding based on multiplexed fragments can be cost-effective and useful for biomonitoring zooplankton in natural communities.

ACK N OWLED G EM ENTS
We thank R. Young and S. Adamowicz for providing the reference sequences for many of the zooplankton species included in this study. We also thank E. Brown

CO N FLI C T O F I NTE R E S T
None declared.

DATA ACCE SS I B I LIT Y
The scripts used for bioinformatics analysis are available in the