Hidden biodiversity in Neotropical streams: DNA barcoding uncovers high endemicity of freshwater macroinvertebrates at small spatial scales

In addition, given recent concerns with the use of the K2P model for species delimitation (e.g., [62], we applied two additional single-locus analyses to confirm lineage diversity: the Bayesian General Mixed Yule Coalescent model (GMYC; [63] and the Automatic Barcode Gap Discovery (ABGD; [64]. The GMYC approach uses branch lengths to determine the tran- sition from intraspecific to interspecific relationships [63]. The ABGD algorithm allows to sort DNA sequences into “hypothetical species” based on the gaps in the distribution of intra- and inter-specific genetic divergence in a given sample [64]. Although the two approaches differ in their properties (i.e., tree branch length vs. distribution gaps), we used them to confirm the patterns of species delimitation. To perform GMYC tree-based analyses, we used the ultrametric trees previously generated with BEAST. GMYC was performed using the single threshold parameter at the GMYC web- server (https://species.h-its.org/gmyc/). ABGD was carried out using the online version of ABGD software [64]; https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html). Default settings were used, however, distance matrices based on K2P distance calculated in MEGA7 were used as input. All analyses were run using a relative barcoding gap width (X value) set to 1.0. Only the recursive results were used because they allowed for different gap thresholds among taxa. To compare patterns of spatial variation in genetic diversity (MOTUs), we quantified the number of shared species among sampling sites. We also estimated Fisher’s alpha index of

diversity and Whittaker’s measure of β- beta diversity. Given that these analyses might be affected by variation in sampling size, we also used rarefaction and extrapolation methods [65,66] as implemented in the R package iNEXT [67]. This method allows for comparisons between sites while controlling for differences in abundance and sampling effort. For these anal- yses, we fit curves for the first three Hill numbers: species richness (q = 0), the exponential of Shannon entropy (“Shannon diversity”, q = 1), and the inverse Simpson concentration (“Simp- son diversity”, q = 2), using individual-based abundance data. Given the limitation in sample size, we did attempt to make statistical inferences regarding differences in diversity among sites. Finally, to assess lineage diversity and the diversity of unique haplotypes across watersheds spanning the Isthmus of Panama, we collected haplotype information (426 haplotypes: Gen- Bank accession numbers KR134410-KR134835) from a previous study in this region [25]. After adding these sequences to our dataset, we generated multiple sequence alignments with MAFFT 7.313 [45] using the L-INS-i algorithm in Geneious V7.03. Then, we trimmed the sequences to the same fragment size and compared the previously reported haplotypes with the ones encountered in our dataset. To exclude redundancies prior to phylogenetic analyses, we applied DAMBE v. 6.4.11 [68] to identify and remove duplicate haplotypes from our data- set. In total, we found a total of 15 duplicate haplotypes between the two datasets. The com-

bined dataset contained 12 sites (8 from Mu´rria et al. [25], and 4 from the present study). One site (Frijolito) was sampled during both studies, but we analyzed them separately to preserve independence between the two studies. We then applied the same rarefaction/extrapolation approach described above to generate rarefaction curves as a function of the number of indi- viduals sampled. Although the collection method was similar (i.e., both studies used kick-nets during a given amount of time), the two studies may not be directly comparable due to differ- ences in sampling effort and the overall study objective. Therefore, we did not attempt to make statistical inferences about the relative abundance of macroinvertebrates among sites or between studies. Instead, we only focused on compiling the total number of unique haplotypes or molecular species that are currently known within each site in this region. In addition, our ultimate goal was to explore the overall molecular endemicity of aquatic macroinvertebrates, rather than providing precise estimates of species diversity in this region. Thus, substantially more research that includes more data, sites, and replication is needed to confirm if current patterns hold across the entire region. PLOS ONE Hidden biodiversity in Neotropical streams PLOS ONE | https://doi.org/10.1371/journal.pone.0231683 August 7, 2020 5 / 13