Detecting fish assemblages with environmental DNA: Does protocol matter? Testing eDNA metabarcoding method robustness

Environmental DNA (eDNA) methods have recently emerged as noninvasive powerful methods to inventory biodiversity in a wide range of ecosystems (reviewed in Tab erlet et al., 2018; Zinger et al., 2020). The recent striking progress in DNA sequencing technologies allowed to shift from single-species detection (eDNA barcoding) (Ficetola et al., 2008) to the simultaneous detection of entire species assemblages (eDNA metabarcoding) (Civade et al., 2016; Received: 19 June 2020 | Revised: 22 October 2020 | Accepted: 23 October 2020 DOI: 10.1002/edn3.158

. A particular interest has been devoted to the development of eDNA metabarcoding studies in aquatic environments (Rees et al., 2014) as water acts as a recipient for DNA, allowing an integrative assessment of aquatic biodiversity Zinger et al., 2020). Among those studies, eDNA fish inventories were particularly effective and often outperform traditional sampling methods (Cantera et al., 2019;Fujii et al., 2019;Pont et al., 2018). Indeed, Cilleros et al., (2019) showed that eDNA fish inventories are not biased by distinct detection probabilities among species or environments whereas traditional capture methods are species selective (e.g., size dependence for gillnets) or inefficient in some habitats (e.g., deep or encumbered areas for electrofishing and seine nets, respectively) (Fujii et al., 2019). The flourishing development of eDNA-based studies led to the development of various eDNA sampling and laboratory protocols, but robust experimental designs allowing between-studies comparisons are still lacking (Bylemans et al., 2018;Dickie et al., 2018;Zinger et al., 2019Zinger et al., , 2020. The methods developed for aquatic eDNA metabarcoding studies generally follow three major steps: (a) eDNA collection and preservation; (b) eDNA extraction; and (c) eDNA detection (Tsuji et al., 2019) with strong variations in each step depending on the studies. In fact, researchers developed various strategies for aquatic eDNA-based studies adapted to field specificities and protocols are rapidly evolving with technology and molecular biology advances.
For instance, the protocol of Valentini et al. (2016) for running waters has widely been used from 2014 to 2016, before being replaced by the protocol of Pont et al., (2018) since 2016 (Figure 1a).
The differences between those protocols lie primarily in the eDNA collection method. Both above cited protocols use filtration with enclosed filters to collect eDNA, this filter type being recognized as more efficient than "open filters" such as disk membranes Spens et al., 2017). Nevertheless, enclosed filter membrane can vary in surface, porosity, and composition. For instance, filters with larger pore size and membrane surface are less sensitive to clogging and allow filtering larger water volumes than those with smaller pore sizes, but might also be less efficient to capture eDNA.
Importantly, such variations can affect the probability of species detection as they are influencing the amounts of captured eDNA (Capo et al., 2019;Li et al., 2018;Takahara et al., 2012;Turner et al., 2014).
In addition, the membrane composition also influences the quantity of eDNA retrieved (Deiner et al., 2018;Hinlo et al., 2017;Majaneva et al., 2018;Shu et al., 2020). The extraction step appears more consensual, with most studies processing eDNA with commercial DNA extraction kit Tsuji et al., 2019;but see Turner et al., 2014 for alternatives). In the same way, most studies agree that during the amplification process, the number of PCR replicates is important to minimize missing detection of taxa that are present (false negatives) (Ficetola et al., 2015) with 3 PCR replicates being the lower limit to get relevant data (Pont et al., 2018;Rees et al., 2014;Valentini et al., 2016) and up to 12 replicates when detection probabilities are low (Ficetola et al., 2015). Similarly, sequencing depth was also claimed to be an important parameter to increase the detection probability of rare taxa (Smith & Peay, 2014). Aquatic eDNA metabarcoding studies hence exhibit considerable variations in the protocols. Considering that those variations can impact species detection probability, it can be expected that different protocols will produce distinct biodiversity measures for the same sampled fauna.
Here, we investigate the robustness of an aquatic eDNA metabarcoding method through the comparison of two biodiversity datasets generated by two different protocols. Those two protocols differ mainly in the eDNA collection step, each protocol using a distinct filter, with different membrane composition, surface, and pore size, leading also to different volumes of filtered water and in the extraction step, with two different extraction kits and different sample preparation prior to the extraction (Figure 1a). Those two protocols were used in many rivers and streams throughout the world including Europe, South America, and Asia (Cantera et al., 2019;Cilleros et al., 2019;Civade et al., 2016;Dufresnes et al., 2019;Lopes et al., 2017;Milhau et al., 2019;Pont et al., 2019;Sasso et al., 2017;Spitzen -van der Sluijs et al., 2020;Vimercati et al., 2020). It is F I G U R E 1 eDNA workflow for the two considered protocols (a) and sampling site locations (b) | 3 COUTANT eT Al. therefore important to test the robustness of eDNA results to such variations in protocols to determine whether the data acquired using distinct protocols deserve to be merged and ultimately if they can, or not, be used jointly for biodiversity and conservation studies. This is of particular importance for the studies replicated over years with temporal changes in the protocols, as well as for data sampled in remote and difficult to access locations.

| eDNA collection
This study was conducted in French Guiana (Latitude 2 to 6°N; Longitude 51 to 54°W). This 80,000 square kilometer territory is part of the Guiana shield and bordered by Brazil to the South and East and by Suriname to the West. Climate is equatorial with annual rainfall ranging from 3,600 mm (North-East) to 2,000 mm (South and West), and most of the territory (90%) is covered by primary forest.
French Guiana also possesses a dense hydrographic network constituted of six major watersheds and several coastal rivers. Altogether, these freshwater bodies constitute a hotspot of freshwater fish diversity and endemism with at least 367 strictly freshwater fish species (Keith et al., 2000). eDNA collection was achieved using two different sampling protocols, protocol 1 (P1) described in Valentini et al. (2016) and protocol 2 (P2) described in Pont et al. (2018). These protocols were developed by SPYGEN and optimized with large volumes of filtered water, high number of PCR replicates and with the amplification of short fragments. P1 and P2 were used to inventory fish fauna from 15 rivers and stream sites in French Guiana ( Figure 1b and Table S1). P1 was used to collect eDNA in 2014 and 2015 by Cilleros et al. (2019). 15 sampled sites were selected for which the bioinformatic pipeline was reprocessed to perform taxonomic assignments based on an updated reference database (containing 24 species more than the one used in Cilleros et al. (2019)). P2 was used to collect eDNA from 2016 to 2019 on the same 15 sites and the taxonomic assignments were performed based on the same updated reference database. Differences between P1 and P2 are provided in Figure 1a, and the general framework of the sampling was as follows.
All the samples were collected in November (dry season), and none of the sites experienced environmental or anthropogenic changes during the study period. The 15 sites (3 large rivers and 12 streams) were sampled using encapsulated filtering cartridges ( Figure 1b and Table S1). In each site, one field replicate was performed by filtering a water volume of c. 50L for P1 (Envirocheck HV, PALL, USA) and c.
34L for P2 (VigiDNA, SPYGEN, France). A peristaltic pump (Vampire sampler, Burlke, Germany) and disposable sterile tubing were used to pump the water through the encapsulated filtering cartridges. The input part of the tube was held a few centimeters below the surface in rapid hydromorphologic units to allow a better homogenization of the DNA in the water column. When filters began to clog, the pump speed was decreased to avoid material damages. To minimize DNA contamination, the operators remained downstream from the filtration either on the boat or on emerging rocks. After filtration, the capsules were filled with a preservation buffer and stored in the dark at room temperature for less than a month before DNA extraction. The three river sites (R 1 to R 3 ) were wider than 30 meters and deeper than 1 meter. The 12 stream sites (S 1 to S 12 ) were less than 15 meters wide and 1 meter deep. The sites belong to different watersheds (Maroni (R 2 , R 3 , S 4 , S 5 ), Oyapock (S 6 ), Sinnamary (S 2 ), Organabo (S 1 ), Mana (S 3, S 8, S 9 , S 11 , S 12 ), Approuague (R 1 ), Comté (S 7 ), and Kourou (S 10 )), and account for various environmental conditions and fish species diversities ( Figure 1b and Table S1). Among sites, turbidity varied from 0.28 to 16 NTU, pH from 4.78 to 7.64, temperature from 22.3 to 31.1°C, conductivity from 15 to 54 µ/S, and O 2 saturation from 45.1% to 111%, encompassing most of the range of freshwater conditions encountered in French Guiana.

| eDNA extraction
Extraction protocols followed different procedures for P1 and P2 for 2 hr, agitated manually for 5 min, and then emptied into three 50 ml tubes. In total, ~120 ml was divided among three tubes that were centrifuged for 15 min at 15,000 g. The supernatant was removed with a sterile pipette, leaving 15 ml of liquid at the bottom of each tube. Subsequently, 33 ml of ethanol and 1.5 ml of 3 M sodium acetate were added to each 50 ml tube. The three tubes were centrifuged at 15,000g for 15 min at 6°C, and the supernatant was discarded. After this step, 360 µl of ATL Buffer of the DNeasy Blood & Tissue Extraction Kit (Qiagen) was added to the first tube, the tube was vortexed, and the supernatant was transferred to the second tube (Tréguier et al., 2014). This operation was repeated for all tubes. The supernatant of the third tube was finally transferred to a 2 ml tube, and the DNA extraction was performed following the manufacturer's instructions. Four negative extraction controls were also realized, amplified, and sequenced in the same way as the field samples to highlight possible contaminations.
With P2, the extraction procedure was achieved according to Pont et al. (2018)'s instructions. For DNA extraction, each filtration capsule was agitated for 15 min on an S50 shaker (cat Ingenieurbüro™) at 800 rpm and then emptied into a 50-ml tube before being centrifuged for 15 min at 15,000 g. The supernatant was removed with a sterile pipette, leaving 15 ml of liquid at the bottom of the tube. Subsequently, 33 ml of ethanol and 1.5 ml of 3 M sodium acetate were added to each 50-ml tube and stored for at least one night at −20°C. The tubes were centrifuged at 15,000 g for 15 min at 6°C, and the supernatants were discarded. After this step, 720 μl of ATL buffer from the DNeasy Blood & Tissue Extraction Kit (Qiagen) was added. The tubes were then vortexed, and the supernatants were transferred to 2-ml tubes containing 20 μl of Proteinase K. The tubes were finally incubated at 56°C for two hours. Afterward, DNA extraction was performed using NucleoSpin ® Soil (MACHEREY-NAGEL GmbH & Co., Düren Germany) starting from step six and following the manufacturer's instructions. The elution was performed by adding 100 µl of SE buffer twice. Four negative extraction controls were also realized, amplified, and sequenced in the same way as the field samples to highlight possible contaminations.

| Amplification and sequencing
DNA amplification and sequencing protocol were the same for P1 and P2. The samples were first tested for inhibition using qPCR following the protocol of Biggs et al. (2015). Then, the samples were diluted 5-fold before the amplification if they were consid- Basel, Switzerland) were also added to the mixture. 12 PCR replicates per field sample were performed. The forward and reverse primer tags were identical within each PCR replicate. The PCR mixture was denatured at 95°C for 10 min, followed by 50 cycles of 30 s at 95°C, 30 s at 55°C, and 1 min at 72°C and a final elongation step at 72°C for 7 min. The amplification step was realized in a dedicated room with negative air pressure and physical separation from the DNA extraction rooms (with positive air pressure). The four negative extraction controls and the three PCR negative controls (also 12 replicates) were sequenced in parallel.

| Bioinformatic processing
The purified PCR products were then pooled in equal volumes to reach a sequencing depth of 400,000 (P1) or 500,000 reads (P2) per sample before library preparation. Five libraries were prepared using the Metafast protocol (https://www.faste ris.com/metafast), a PCRfree library preparation, at Fasteris facilities (Geneva, Switzerland).
Sequencing were performed using an Illumina HiSeq2500 For both protocols, the sequence reads were analyzed using the functions of the OBITools package (http://metab arcod ing.org/ obitools) following the protocol described in Valentini et al. (2016).
The ecotag function was used for the taxonomic assignment of mo- occurring with a frequency bellow 0.0003 per library sample were considered as tag-jumps and discarded (Schnell et al., 2015). These thresholds were empirically determined to clear all reads from the extraction and PCR negative controls included in our global data production procedure as suggested in De Barba et al. (2014).

| Data analyses
According to the river continuum concept, the biodiversity markedly differs along the upstream-downstream gradient (Vannote et al., 1980). This has also been demonstrated in tropical regions and particularly in Guianese streams (Cilleros et al., 2017). As our sites represented a large range of watercourse types and sizes, we ordered the sites using the Strahler stream order. The Strahler stream order defines stream size based on their hierarchy of tributaries, with increasing orders from the source to the estuary (Strahler, 1957). Among the different taxa obtained, only taxa assigned to the species level were considered. Species by site matrices with presence/absence data were created separately for the sites sampled with P1 (P1-samples) and P2 (P2-samples) (Table S2 and Table S3). The differences in species richness and species occurrence (the number of sites where each species was detected) between the P1-and P2-samples were tested with a Mann-Whitney U-tests for nonparametric and unpaired samples (both protocol having been conducted at different years). The dissimilarity (β-diversity) was calculated between sites within each protocol for each habitat separately (for each protocol: 66 sample pairs for the streams, 3 sample pairs for the rivers), and between P1-P2 sample pairs (15 pairs) of each site with the Jaccard dissimilarity index (β jac ) using the betapart package (Baselga & Orme, 2012) and the functions beta.pair and beta.temp. The two components of β jac , nestedness (β nes ) and turnover (β turn ) as well as their individual contribution percentage (P turn and P nes ) were then calculated for the dissimilarity values between P1 and P2 site pairs as follows: P turn / β jac × 100 and P nes /β jac × 100 (Villéger et al., 2014). This quantifies the extent to which the species composition (turnover) and richness (nestedness) contribute to differentiate the communities of each sample pair for each site (Baselga et al., 2018).
Protocol comparisons were led by considering β jac dissimilarity values at two different levels ( Figure 2): (i) between sites within each protocol and, (ii) between P1 and P2 sample pairs of each site. The dissimilarity values of the first level (i) were compared between both protocols using Mann-Whitney U-tests ( Figure 2A) to test whether the two protocols present replicability differences in both habitats (stream and river considered separately). Likewise, the dissimilarity values of the first level (i) were also compared between habitats using Kruskal-Wallis tests ( Figure 2B) to determine whether both protocols have the same performance in stream and river (protocols considered separately). Finally, the dissimilarity values of the second level (ii) were compared between both habitats using Kruskal-Wallis tests ( Figure 2C) to determine whether sample dissimilarity between P1 and P2 is dependent on habitat type. In addition, a unilateral one sample Z test was performed to test whether turnover contribution to β jac dissimilarity was higher than 50% (habitat type considered separately), and therefore whether turnover and nestedness relative contribution to P1-P2 samples dissimilarity are dependent on habitat type.
As the β-diversity and the α-diversity are interrelated, the β jac is influenced by the number of species present in each pair of samples. To address this issue, a null model using the Raup crick metric (β rc ) was performed to determine whether the pairs of samples were more similar than expected under stochastic effects. This metric expresses the dissimilarity between two communities relative to the null expectation of randomly assembled communities by estimating the probability that two randomly drawn communities have the same number or more species in common than the two observed communities. The β rc metric was calculated between sites within each protocol for each habitat separately (for each protocol: 66 sample pairs for the streams, 3 sample pairs for the rivers), and between P1 and P2 sample pairs of each site following Chase et al.. (2011). For each site, 9,999 mock communities were simulated for each sample pair by randomly selecting species from the regional species pool (all the species detected in the samples). Each simulated community contained the same species richness to that of the observed community, and the presence probabilities of each species were weighted by its among-site occupancy. The dissimilarity was calculated for each pair of simulated communities and the observed β-diversity was then compared to the null distribution, e.g., β rc = β obs -E(β null ) where β obs is the observed dissimilarity of the observed sample pair communities and E(β null ) is the mean dissimilarity of the simulated sample pair communities. The resulting value represents a β-diversity metric ranging from −1 to +1. A negative β rc value stands for communities more similar than expected by chance while a positive value refers to communities less similar than expected by chance. A null value indicates that community assembly is highly stochastic.
The β rc dissimilarity values were considered at the same first level (i) and compared with the same statistical tests as detailed above for the β jac analysis (Figure 2). The second level (ii) of comparison was investigated using a permutational analysis of variance (PERMANOVA; Anderson, 2001) to test for differences in species composition after checking for homogeneous dispersion among P1-and P2-samples with a permutation dispersion (Anderson et al., 2006). Nonmetric multidimensional scaling (NMDS) were then performed on the Jaccard and Raup-Crick indices to visualize the stochastic and nonstochastic effects explaining sample pair dissimilarities. All the analyses were performed using R software version 3.6.1 (2019-07-05) (R Core Team, 2019).

| RE SULTS
A total number of 2,251,622 and 2,430,314 reads were retained after bioinformatic screening for the sites sampled with P1 and P2, respectively. They, respectively, represented 48.09% and 51.91% of the 4,681,936 total reads retained. After the bioinformatic filtering step, reads were detected in all our 360 PCR replicates while no reads were found in the PCR and extraction controls. Among all the samples, we detected 160 species out of the 255 species referenced in our database. We obtained 779 occurrences in the 15 sites, and all the species detected were consistent with their known distribution.

| Protocol effect on species richness and occurrence
Among the 15 sites, we detected 139 species with P1 and 148 species with P2. The species richness detected tended to increase with the watercourse size of the sites, regardless of the protocol (Figure 3a), F I G U R E 2 The different levels of dissimilarity comparison. (i) Between site dissimilarities obtained with P1 were compared to that of P2, streams and rivers were tested separately (A); between stream site dissimilarities were compared to between river site dissimilarities, P1 and P2 were tested separately (B). (ii) P1 and P2 sample pair dissimilarity of streams were compared to that of rivers (C). Circles with the same colour refer to the same sites sampled with P1 or P2, dashed lines represent the dissimilarities considered, (i) / (ii) indicates the level at which dissimilarity values were considered, arrows show the different comparisons (A, B and C; see Methods) emphasizing the difference in species richness between streams and rivers. The low species richness detected in the P1-sample of the R 1 site is suspicious and may reveal an issue during the sampling or laboratory treatment (Figure 3a). Nevertheless, excluding site R1 from our data did not changed our conclusions (results not shown). Among the 15 sites, the number of detected species with P1 ranged from 4 to 80 (median = 22) and from 3 to 68 (median = 27) with P2. The number of species detected with P2 was higher than with P1 in 9 out of the 15 sites, but the average species richness obtained among sites did not significantly differ between protocols (Mann-Whitney U rank test, U = 86.5, p = 0.289, n = 15) ( Figure 3b).
Among the 779 occurrences obtained across all the sites, 341 and 438 occurrences were detected in the P1-samples and the P2samples, respectively. Moreover, among the 160 species detected, 12 species were only detected in the samples collected with P1 while 21 species were only found with P2. Among the 12 species recovered only with P1, 10 occurred in a single site and 2 species were present in 2 sites. Among the 21 species found only with P2, 16 occurred in a single site, 4 in two or three sites, and one species was detected in 4 sites. Regarding the 127 species detected with both protocols, most of the species occurred in 1 or 2 sites, representing 69.29% of the P1-and 59.06% of the P2-species, respectively ( Figure 4a). A few of the shared species occurred in more than half of the sites with higher occurrences detected with P2, representing 7.09% and 9.45% of the P1-and P2-species, respectively. The species obtained with both protocols represented 79% of the total fauna detected. The species occurrence ranged from 1 to 10 (median = 1) in the P1-samples while it ranged from 1 to 14 (median = 2) for the P2-samples. Overall, there was a significant difference in the average number of species occurrence between protocols, with P2 detecting from 0 to 4 more occurrences than P1 (Mann-Whitney U rank test, U = 15,056, p = 0.004, n = 160) (Figure 4b).

| Protocol effect on community composition
Considering the first level (ii) of comparison (between site dissimilarities within each protocol), the β jac and the β rc dissimilarity of the within-P1-samples did not significantly differ to that of the within-P2-samples neither for the rivers nor for the streams (Figure 2A and Table 1). This indicates that the degree of species dissimilarity among sites is not significantly different between protocols in both habitats.
Moreover, the β jac dissimilarity of the within-P1-stream samples did not significantly differ to that of the within-P1-river samples ( Figure 2B, Kruskal-Wallis rank-sum test, χ 2 = 3.28, df = 1, p = 0.07, n = 69). Likewise, the β jac dissimilarity of the within-P2-stream samples did not significantly differ to that of the within-P2-river samples ( Figure 2B, Kruskal-Wallis rank-sum test, χ 2 = 0.15, df = 1, p = 0.70, n = 69). This indicates that both protocols perform in the same way in rivers and in streams when considering both nestedness and turnover components of β-diversity. However, the β rc dissimilarity values of the within-P1-and the within-P2-stream samples were both F I G U R E 3 Site species richness obtained with protocol 1 (P1) and protocol 2 (P2) (a) and averaged across sites for each protocol (b). The Strahler stream orders are indicated below the bars of panel (a), and the sites are ordered by decreasing stream order (indicating the size of the river/stream, with larger ones having higher Strahler orders. Bright orange and dark grey bars indicate the sites sampled with P1 and P2, respectively. Species richness of P1 and P2 samples were compared using a Mann-Whitney U-test, NS = non significant (p > .05)
Considering the second level (ii) of comparison, the dissimilarity difference between P1 and P2 stream sample pairs of each site and P1 and P2 river sample pairs of each site, the β jac dissimilarity ranged from 0.28 to 0.93 (median = 0.29) in the three river sites and from 0.25 to 0.91 (median = 0.52) in the twelve stream sites. There was no significant difference between the dissimilarity of the P1-P2 river sample pairs and the P1-P2 stream sample pairs ( Figure 2C, Kruskal-Wallis rank-sum test, χ 2 = 0.08, df = 1, p = 0.87, n = 15) ( Figure 5a). In addition, nestedness contribution to sample pair dissimilarity ranged from 13.24% to 100% (median = 41.78) for the river sample pairs while it ranged from 5.38% to 100% (median = 46.36) for the stream sample pairs. For the sites R 1 , S 1 , S 5 , and S 6 , nestedness contributed to 100% of the total sample pair dissimilarity (Figure 5a,b). Turnover contribution to the dissimilarity between sample pairs for each site ranged from 0% to 86.76% (median = 58.22) for the river sample pairs while it ranged from 0% to 94.62% (median = 53.64) for the stream sample pairs. However, the turnover contribution to site dissimilarity was not significantly greater than 50%, thus implying no significant difference between the nestedness and turnover contribution to β jac dissimilarity (one sample Z test, z = 3.60, p = 1, n = 15).
Moreover, the turnover contribution as well as the nestedness contribution did not significantly differ between river and stream sam- two distinct clusters. Moreover, the R 1 samples were the most distant samples belonging to the same site due to the large difference in species richness obtained between P1 and P2. The sample pairs of the S 6 and S 9 sites were also more distant from each other than the other sample pairs (Figure 6a). This is consistent with the previous results demonstrating that the S 6 sample pair presented the highest dissimilarity due to nestedness and that the S 9 samples exhibited the highest β jac dissimilarity among stream sites (Figure 5a,b). The use of the null model allowed correcting for the difference in α-diversity and therefore only account for compositional differences in the communities between sample pairs. The resulting β rc index indicated that in all sites, the P1-and P2-samples presented homogeneous dispersion While controlling for species richness differences in the β rc NMDS, the river samples remain more distant from each other than the stream sites. Specifically, the closeness between R 1 , S 6 , and S 9 sample pairs in the β rc NMDS indicate that their strong dissimilarity identified by the Jaccard's index was mainly due to differences in the number of detected species, while species turnover remained low (Figure 6b).

| D ISCUSS I ON AND CON CLUS I ON
The variety of protocols developed for aquatic eDNA-based metabarcoding studies represents a hindrance to study comparisons (Minamoto et al., 2016;Shu et al., 2020). Moreover, because F I G U R E 5 β jac dissimilarity (a) and nestedness and turnover contribution to β jac dissimilarity (b) of sample pairs for each site. White and striped bars represent nestedness and turnover, respectively. Sample pairs refer to samples collected with P1 and P2 in the same site (a) (b)

F I G U R E 6
Non-metric multi-dimensional scaling (NMDS) ordinations based on β jac (a) and β RC (b) indices. The segments connect each pair of samples from the same site collected with P1 and P2. River sites and stream sites are clustered with a grey area. Communities that are closer together, using modified Raup-Crick dissimilarities, are more deviant from the null expectation, whereas communities that are farther apart are less deviant from the null expectation protocols are continuously improved and optimized with scientific advances, heterogeneous datasets are produced, some of which may become deprecated, for either technical or financial reasons (Deiner et al., 2015). Understanding the extent to which different protocols produce different outputs by comparing community species richness and composition is therefore a key question to determine whether merging data obtained with distinct protocols is possible, or whether the biodiversity estimates from distinct protocols require to be considered separately. The two protocols compared in this study followed the same general optimized workflow but presented marked differences during the data collection and the extraction step.
Overall, the two protocols produced similar measures of richness in species. Although P2 allowed, on average, to detect more species per site than P1, such differences remained low and nonsignificant when considering the sites altogether. From an ecological point of view, the observed richness obtained with the two different protocols remains consistent with the expected richness of fish fauna across Guianese rivers (Le Bail et al., 2012). Moreover, the increasing richness pattern along the upstream-downstream continuum revealed by both protocols supports the river continuum concept, which describes changes in species assemblages along the upstream-downstream gradient in response to variability in the environment and a gradual increase of the stream size offering more habitats (Vannote et al., 1980). This richness pattern also echoes Cilleros et al. (2017) findings that demonstrated an increase in local species richness with stream size from the upstream to the downstream of the same Guianese rivers using traditional sampling data.
The weak richness of the R 1 sample collected with P1 might thus reveal an issue during the sampling or the laboratory processes as it is the only river sample presenting a lower richness than the stream samples.
The known patterns of species occurrences were conserved in the samples regardless of the protocol, even if P1 was slightly less efficient at detecting the most frequent species known to be widespread in Guianese streams and rivers such as Copella carsevennensis, Hoplias malabaricus, or Gymnotus coropinae (Planquette et al., 1996).
Species inhabiting most of the small streams such as Nannacara aureocephalus and Callichthys callichthys (Cilleros et al., 2017;Keith et al., 2000) were also rightly retrieved among the most occurring species in both P1-and P2-stream samples. Overall, the occurrence patterns obtained with both protocols were congruent with the known proportion of fish with large and small distribution encountered in Guianese rivers. Indeed, we recovered around two third of fish species with restricted distribution (occurring in one or two sites) and one third of fishes with wider geographical distribution, paralleling previous studies on the spatial distribution and species composition of the Guianese fish fauna (Cilleros et al., 2016;Le Bail et al., 2012). More generally, our occurrence gradients are in line with the hyper-diversified tropical regions recognized to present a majority of species with a restricted spatial distribution and only a limited number of widely distributed species (Colwell et al., 2004).
P1 and P2 also performed similarly on retrieving patterns of beta-diversity. Both β jac and β rc dissimilarities of within-P1-samples were not significantly different to that of within-P2-samples in both streams and rivers indicating that the protocols allowed for similar replicability in both habitats. In addition, the protocols provided similar species inventories for each site regardless of habitat type, as no significant compositional difference among samples collected with P1 and P2 was found. Moreover, the absence of a significant difference in β jac dissimilarity between river and stream samples within each protocol suggested that both protocols performed as well in streams as in rivers. This was also supported by the absence of significant difference in turnover contribution to β jac dissimilarity between P1 and P2 stream sample pairs and P1 and P2 river sample pairs of each site.
However, the smaller β rc dissimilarity values in within-P1-and within-P2-stream sample pairs than in within-P1-and within-P2river sample pairs, respectively, indicated that river samples presented more differences due to species composition than stream samples. This was also emphasized by the larger size of the cluster formed by the river samples on the β rc NMDs than on the β jac NMDS.
This demonstrates that when considering only species turnover and fixing for the nestedness effect, P1 and P2 protocols allow the same sample replicability but they both perform better in streams than in rivers. This is in line with Cantera et al. (2019) (12), which might also contribute to the robustness of the method. Indeed, Ficetola et al. (2015) emphasized the importance of adjusting the number of PCR replicates depending on the features of the study system and suggested to realize at least 8 PCR replicates when dealing with taxa with low detection probabilities (e.g., rare species). Our large number of PCR replicates was therefore adapted to tropical environments containing rich species assemblages with a majority of rare species.
The general call for developing a standard approach for aquatic eDNA-based metabarcoding studies is challenging because of the specificity of each environment (Deiner et al., 2018). Still, while it has been proven that the protocols and the material used have an impact on the total quantity of DNA recovered (Capo et al., 2019;Thomas et al., 2018), it does not necessarily implies that the different datasets generated by distinct protocols are strongly affected by those protocol variations (Djurhuus et al., 2017). In fact, protocol variations may not produce heterogeneous datasets if the achieved sampling effort provides a valuable representation of the communities. Despite the evident need to search for a unique standardized approach offering the best biological results, testing the robustness of existing methods remains a way to potentially rescue biological data collected using deprecated protocols. This is of particular importance to consider time series data, or biological inventories gathered in remote locations difficult to access.

DATA AVA I L A B I L I T Y S TAT E M E N T
The site locations (Table S1) and the species by sites matrices (