The larvae used in this experiment come from another study using DNA barcoding for species identification. Ichthyoplankton sampling was carried out in the Mogi-Guaçu river, an important tributary of the left bank of the Grande River, Upper Paraná River Basin. Samples were taken from two sites: 1: 21°35.441’S, 47°57.244’W (Dec/16 and Jan/17) and site 2: 21°30.168’S, 48°2.534’W (Dec/15, Nov/16, Dec/16, and Jan/17) (Table 1), located in the northeast region of the state of São Paulo, Brazil. The sites are about 15 km apart from each other in a straight line.

Table 1.Download as CSV 

Number of larvae in each sample and their respective sites and month of capture.

Sample number	Number of larvae	Site	Month
S1	36	1	Jan/17
S2	42	1	Dec/16
S3	26	2	Dec/16
S4	38	2	Jan/17
S5	44	2	Nov/16
S6	44	2	Dec/15
Total	230

Collecting took place with a conical-cylindrical net with a mesh size of 0.5 mm. Samples were immediately fixed in 96% ethanol. In the laboratory, eggs and larvae were separated from other materials present in the samples (e.g., leaves, sediment, and sticks) under a stereomicroscope and kept in 96% ethanol until the molecular procedure.

In this study, we used 230 larvae distributed randomly across six arrays, with the number of organisms per array ranging from 36 to 44. However, each sample contained larvae of only one site and month (Table 1). Each larva was cut into halves; one half was subject to gDNA extraction and subsequent Sanger sequencing, while the other half was pooled and homogenized with other larvae from the same sample for combined DNA extraction and subsequent metabarcoding using next-generation sequencing (NGS) (Figure 1).

Figure 1.

Experimental setup of this study for Sanger and metabarcoding analysis of ichthyoplankton.

DNA extraction for Sanger sequencing was made following the protocol proposed by Ivanova et al. (2006). Cytochrome C Oxidase Subunit I (COI) primers FishF1 and FishR1 (Ward et al. 2005) were selected because these have a good species recovery in the hydrographic basin where the study occurred. PCR conditions were: 1 min at 96 °C followed by 35 cycles of 10 s at 96 °C, 5 s at 50 °C and 4 min at 60 °C. The PCR product was sequenced using a Big Dye 3.1 Terminator kit (Applied Biosystems) in an ABI Prism 3130 (Applied Biosystems), generating two complementary sequences of 654 bp.

DNeasy Blood and Tissue Kit (Qiagen) was used for DNA extraction from bulk samples following the manufacturer’s protocol. Two PCR runs were performed: the first to amplify the COI region, using the same primer pair used for the Sanger sequencing coupled with Illumina MiSeq adapter, and the second to attach a multiplex identifier (MID).

In the first PCR, each reaction contained, 2 μl of DNA template, 2.5 μl buffer, 1 μl MgCl₂ (50mM), 1 μl dNTP (2 mM), 0.5 μl forward adapter+FishF1 primer (5 mM), 0.5 μL reverse adapter+FishR1 primer (5 mM), 0.4 μl Taq Polymerase (5 U/µl) (Phoneutria), and 17.1 μl of molecular biology grade water in a final volume of 25 μl. PCR conditions were 95 °C for 3 min, 35 cycles of 95 °C for 30 s, 52 °C for 45 s, and 68 °C for 60 s; and a final extension at 68 °C for 10 min. We purified PCR products through NucleoSpin Gel and PCR Clean-up from Macherey-Nagel. In the second PCR, the library construction was made, where amplicons were dual indexed with MID’s by Macrogen, using the Nextera XT Index Kit and following the instructions present in 16S Metagenomic Sequencing Preparation. Amplicons were sequenced on the Illumina MiSeq platform, generating paired-end reads of 301 nucleotides in length.

For sequences generated with the Sanger method we used the Geneious Pro v4.8.5 software (Kearse et al. 2012), where forward and reversed files were combined to assemble a consensus sequence. All these consensus sequences were then exported as fasta files to be matched against the same database comprising the Illumina reads.

Metabarcoding data were analyzed following a bioinformatic protocol described in the Suppl. material 1. As the overlapping is not possible due to the length of PCR product (654 bp), only the forward file was used (Cruaud et al. 2017). This file was selected by checking the sequencing quality with FASTQC (Andrews 2010) (Suppl. material 2). Next, a quality trimming with Trimmomatic (Bolger et al. 2014) was performed, and we filtered out reads with a minimum quality of Q20 in at least the first 250 nt. We also cut the sequences up to the length of 250 nucleotides since the quality of final region was low (Suppl. material 2). For this trimming we used the Fastx-Trimmer software (http://hannonlab.cshl.edu/fastx_toolkit/).

After this step, we used the OBITools toolkit version 1.01 22 (http://metabarcoding.org/obitools; Boyer et al. 2016) for DNA metabarcoding analysis. We first used the “obiuniq” command to dereplicate the sequences, to reduce the size file and computation time. In this step, all of the identical DNA sequences were grouped into a single sequence, scoring their abundance for comparison with the Sanger data. Subsequently, we used the “obigrep” command to remove sequences with less than 50 reads.

We constructed a custom database by downloading sequences from Project – FUPR Fishes from Upper Parana River, Brazil, present in the BOLD System (Ratnasingham and Hebert 2007), since this project was executed in the same river basin as the current study. Also, we downloaded sequences from RefSeq (https://www.ncbi.nlm.nih.gov/refseq/), using keywords “fish”, “fishes”, and “COI”, and then filtered the results only for “bony fishes”. For the first step of identification we used the BLASTN software (Altschul et al. 1990) to indicate which sequences had the greatest similarity to those obtained via Sanger and metabarcoding.

As a final filtering step, all sequences generated in the present study (NGS + Sanger) were aligned with those sequences present in a custom database that had similarity greater than 95% with one of the sequencing methods, using the Muscle (Edgar 2004) aligner implemented in Geneious 4.8.5 (Kearse et al. 2012). Next, we calculated the genetic distances among and within the species, using the Kimura-2-Parameter (K2P) distance model (Kimura 1980). To provide a graphic representation of the patterning of divergence among species, a neighbor-joining (NJ) dendrogram of K2P distances was created using the MEGA v 7.0 software (Kumar et al. 2016) for each sample. Those sequences (NGS and Sanger) that presented divergence greater than 2% in the K2P analysis compared to the sequences downloaded from the BOLD System database were removed from the analysis of abundance (Pereira et al. 2013).

To validate the efficacy of metabarcoding, statistical analyses were conducted to evaluate whether this method has a quali-quantitative resolution. To achieve our qualitative aim, from the genetic distance (K2P) and neighbor-joining (NJ) dendrogram, we estimated richness (S = number of species) per sample in each method. Then, the richness between methods was compared to estimate the rate of species detection by metabarcoding, using the Sanger data as a reference.

To achieve our quantitative analysis, as the exact number of larvae per sample was known, the relative frequency per species (or taxon) was estimated for the Sanger data. This frequency compared with the value of the relative frequency (number of reads assigned to a taxon or species in relation to the total) obtained metabarcoding. This comparison was made per sample individually (e.g., Sanger sample 1 × NGS 1 sample). After this, we applied a Pearson correlation to test if the frequencies of Sanger and NGS showed a significant positive correlation and a permutational analysis of variance (Permanova) to evaluate significant differences between relative abundances of Sanger and NGS data.

Of the 230 larvae used in this experiment (Suppl. material deposited in: https://doi.org/10.6084/m9.figshare.6726956) 226 were retained for the analysis with Sanger sequencing. Four larvae were removed, due to low quality in Sanger sequencing and genetic divergence above 2% in the alignment with those downloaded sequences. Results from the Sanger sequencing showed a richness of 29 taxa, belonging to 12 families and three orders, varying from six species in sample S2 to 13 species in sample S1 (Table 2). Only three taxa were assigned to family or subfamily level.

NGS yielded 2,511,656 reads in the eight samples (Suppl. material deposited in: https://doi.org/10.6084/m9.figshare.6726956) and, after the bioinformatic analysis, 633,224 reads were kept. In total, 28 taxa were identified with metabarcoding NGS out of the 29 taxa identified by the Sanger sequencing. Only Astyanax schubarti was not identified by metabarcoding in any of the samples (Table 2).

Table 2.Download as CSV 

Species composition and relative abundance of analyzed samples. S= Sanger sequences; NGS= Next generation sequencing metabarcoding; * false negative.

Family	Species	S1		S2		S3		S4		S5		S6
		S	NGS	S	NGS	S	NGS	S	NGS	S	NGS	S	NGS
Characiformes
Anostomidae	Leporellus vittatus							5.6	10.4
	Leporinus friderici	13.9	23.2					2.8	1.5	9.3	11.6	2.3	1.2
	Leporinus paranensis											2.3	3.7
	Leporinus striatus	5.6	2.5
	Megaleporinus obtusidens			7.3	*	3.8	3.3			9.3	0.8	2.3	0.1
	Megaleporinus piavussu							2.8	0.3
	Schizodon nasutus			2.4	2.1
Characidae	Astyanax lacustris	11.1	0.1					5.6	0.3	4.7	*
	Astyanax schubarti	2.8	*					5.6	*
	Cheirodontinae	11.1	4.3	2.4	0.05	3.8	0.6	2.8	2.0	4.7	1.6
	Cheirodontinae 1					3.8	6.2
	Hyphessobrycon eques							5.6	2.0
Curimatidae	Steindachnerina insculpta	2.8	0.6
Erythrinidae	Hoplias malabaricus									2.3	2.2
Parodontidae	Apareiodon affinis			2.4	0.1
Prochilodontidae	Prochilodus lineatus							2.8	15.2	4.7	9.2
Gymnotiformes
Sternopigydae	Eigenmannia trilineata											2.3	5.0
Siluriformes
Doradidae	Rhinodoras dorbignyi	2.8	1.4
Heptapteridae	Heptapteridae	5.6	2.2							2.3	1.9
	Iheringichthys labrosus											2.3	4.7
	Pimelodella meeki	5.6	2.4					8.3	2.1
Loricariidae	Rhinelepis aspera	2.8	16.9
Pimelodidae	Pimelodus maculatus	27.8	36.0	70.7	95.0	69.2	58.1	52.8	61.3	32.6	49.0	86.4	85.2
	Pimelodus microstoma	2.8	6.3			11.5	21.1
	Pinirampus pirinampu					3.8	3.9
	Pseudoplatystoma corruscans			14.6	2.7	3.8	6.9			16.3	8.6	2.3	0.1
	Pseudoplatystoma reticulatum									4.7	1.3
	Sorubim lima							2.8	3.0	9.3	13.8
Pseudopimelodidae	Pseudopimelodus mangurus	5.6	4.2					2.8	1.8

No false positives were observed in this study (Table 2); all the taxa retrieved by metabarcoding were also found in Sanger sequencing data. However, in four out of the six samples, one taxon was not detected by NGS: S1 and S4 – Astyanax schubarti, S2 – Megaleporinus obtusidens, S5 – Astyanax lacustris. In the other two samples, all taxa detected by Sanger sequencing were also detected by NGS (S3 and S6). The taxa detection rate was higher than 83% in all samples, with 100% of detection in four samples and higher than 90% in three others (mean, 95%).

The correlation between Sanger and NGS methods showed few differences in all samples. We detected a considerable difference between the two methods in only some cases: Rhinelepis aspera, Cheirodontinae sp.1, and Astyanax lacustris, in sample 1 and Prochilodus lineatus and Astyanax lacustris in sample 4. Nevertheless, the samples showed non-significant differences between the relative abundances of the two methods in PERMANOVA analysis (Pseudo-F = 0.39293; P(perm) = 0.838), denoting that, despite some deviations in the values, the quantification of organisms with the NGS method can be considered reliable here (Table 3; Figure 2). Also, in the same way, all samples showed higher values of Pearson’s correlation in relative frequency between Sanger and NGS methods, varying from 0.776 (sample S1) to 0.998 (sample S6) (Table 3; Figure 2).

Figure 2.

Graphic representation of species abundance proportion (%) between Sanger (x-axis, left side) and NGS (x-axis, right side) for analyzed samples. A: sample 1, B: sample 2, C: sample 3, D: sample 4, E: sample 5, F: sample 6.* = false negative.

Table 3.Download as CSV 

Species richness per sample and method of sequencing, detection rate, and Pearson correlation (r) between Sanger and NGS methods.

Sample	Richness		Detection rate NGS/Sanger	Pearson correlation
	Sanger	NGS
S1	13	12	92.3	0.776
S2	6	5	83.3	0.987
S3	7	7	100.0	0.976
S4	12	11	91.7	0.955
S5	11	10	90.9	0.914
S6	7	7	100.0	0.998

Despite some failures in species detection, we can consider that the use of NGS is functional for COI-based ichthyoplankton species detection, with an average detection rate higher than 95%, which is similar to previous findings that compared Sanger and NGS (Cruaud et al. 2017). Also, when compared to morphology-based methods of egg and larvae identification, detection rates were even more satisfactory, as more than 50% of the larvae caught in natural populations cannot be identified due to the absence of diagnostic characters (newly hatched larvae) or damage due to fixation process (Bialetzki et al. 2005; da Silva et al. 2015; Silva et al. 2017).

With our present approach, most individuals were identified to the species level, demonstrating that even using only the forward reads, the information was sufficient for species identification. Only three taxa were assigned to family or subfamily level only, probably because they are cryptic species, not described or with sequences not deposited in any database. This observation reinforces the idea that more sequences of fish species must be deposited in databases, especially in river basins with rich fish fauna as the Upper Paraná River. Another positive aspect is that, compared with DNA barcoding (Arroyave and Stiassny 2014, Becker et al. 2015, Frantine-Silva et al. 2015), our method can trigger significant advances in the quality of ecological studies in fish, as observed in some recent studies with metabarcoding (Guillerault et al. 2017, Maggia et al. 2017, Mariac et al. 2018, Barbato et al. 2019) and, mainly, to ensure significant reduction in labor time and cost per specimen reaching 80% of reduction (Shokralla et al. 2015).

Despite some differences in obtaining DNA samples, both environmental DNA (eDNA) and DNA metabarcoding from bulk samples are subject to false negatives and false positives (Ficetola Gentile et al. 2015, 2016; Divoll et al. 2018). However, in the present study, only false negatives were observed. False negatives can lead to losses in mitigating measures of threatened species, due to the non-detection of target species, which might lead to underestimation of populations and non-implementation of mitigation measures, causing injury to these species (Ficetola Gentile et al. 2015, 2016; Divoll et al. 2018). False negatives are mainly those related to inhibition issues, primer bias, PCR or sequencing errors, and bioinformatics analysis (Ficetola Gentile et al. 2015, 2016; Divoll et al. 2018).

Quali- and quantitative studies involving metabarcoding and eDNA have increased significantly in recent years (Hänfling et al. 2016; Valentini et al. 2016; Evans and Lamberti 2017; Yamamoto et al. 2017; Nakagawa et al. 2018; Mariac et al. 2018). However, the reliable quantification of organisms remains a problem (Evans et al. 2016; Kelly et al. 2017; Hansen et al. 2018; Mariac et al. 2018). eDNA studies have shown a positive correlation between species-specific DNA concentrations in biomass and abundance in eDNA (Takahara et al. 2012; Maruyama et al. 2014; Doi et al. 2015; Evans et al. 2016). In this way, by focusing on developing a low-cost yet efficient protocol, here we show that PCR-based methods can be used to determine relative abundance in studies of eggs and larvae.