Unveiling biogeographical patterns of the ichthyofauna in the Tuichi basin, a biodiversity hotspot in the Bolivian Amazon, using environmental DNA

incompleteness of the COI database for Amazonian fishes, this barcode is currently the one with the best coverage in the number of species. Blastn results were parsed using MEGAN software [44] in order to perform affiliation with a minimum score value of 150, a minimum percent identity of 97% and the naive LCA algo- rithm at minimum percent coverage of 50. The three types of negative controls (field, extrac- tion and PCR) were used to identify and remove contaminating reads in sample libraries with the R package Microdecon [45], which uses the data in the 42 negative controls performed throughout the entire project to identify and exclude the reads in each library that correspond to contamination. The Maximum Sensitivity plus Specificity (MaxSSS) threshold value [46] was computed using the two positive controls (species composition known) and the R package ROCR 1.0–7 [47]. This minimum frequency at which a taxon is retained was applied in order to maximize the sum of True Positive Rate and True Negative Rate. To assess whether the sequencing effort was sufficient, the rarefaction curves for each library were performed with the R Vegan 2.5.6 package [48]. Taxonomic diversity and geographic distribution. Taxa richness (alpha diversity) and com- munity composition (beta diversity) were computed with R packages phyloseq 1.28 [49] and heatmap3 1.1.6 [50]. To detect possible biogeographic entities, a hierarchical cluster analysis (HCA) was per- formed using the hclust function of the vegan 2.5.6 R package [48] to build a dendrogram clus-

tering the sites according to their level of faunistic identity using Jaccard’s binary distance and a Ward2 aggregation method. A bootstrap test with a number of 1000 replications was per- formed to evaluate the robustness of each node of the dendrogram using the package pvclust 2.2.0 [51]. A precise identification of a biogeographic structure with eDNA samples collected in the aquatic environment is only possible if its diffusion is spatially limited. We have there- fore estimated, among our samples, the minimum distance at which species present at one site are not detected at the downstream site. In order to assess the relevance of our sampling effort we computed the second order jac- knife (Jack2) estimator [52] using the package Vegan with the functions specnumber and spec- pool and calculated the percentage of coverage as: observed number of species / Jack2 estimator) calculated for the entire study area and for each biogeographic entities defined by hierarchical cluster analysis. It should be noted that the Jack2 estimator (derived from the number of singletons and doubletons observed) only partially represents the actual species richness of the study area due to the incompleteness of the current COI database which con- tains only 42.6% of the known species in this Amazon watershed. Moreover, the species rich- ness revealed in this study is also partial because our sampling was only carried out during the flooding period and therefore lacks species whose presence would be exclusively seasonal and outside this period.

Environmental factors correlated to the biogeographic structure. The link between environ- mental factors and the biogeographic structure of fish species was measured for four categories of variables noted at each site: 1) the geomorphological environment of the habitat, including the lake or river ecotype, the size of the lakes and their age, as well as the area of their hydro- graphic basin and the length of the hydrographic network that drains into them, 2) the geo- graphical location, including latitude (LAT), longitude (LONG) and altitude (ALT), 3) the substrate, corresponding to the nature of the sub-soil and soil at the sampling point, according to five geological types (GEOL) and six pedological types (PEDO) and, 4) hydro-ecoregions (HER), which integrate some of the above variables and define 7 different categories for the sites sampled in this study (Fig 1). Lake sizes, hydrological networks length and basin area val- ues were extracted for each site using Hydro Atlas [53] under QGIS version 3.4.12 [54] from GIS layers obtained at the World Data Bank (https://datacatalog.worldbank.org/dataset/ PLOS ONE Fish diversity and distribution in the Tuichi basin PLOS ONE | https://doi.org/10.1371/journal.pone.0262357 January 4, 2022 7 / 25