DNA Barcode Reference Library and Undetected Diversity of Fish Species in the Yuanjiang River, China

1. Introduction

The Red River is an internationally important river that flows through China and Vietnam. The upper reaches of the Red River, the Yuanjiang River, is an essential component of the Mountains of Southwest China biodiversity hotspot, which harbors rich fish diversity and endemic fish species due to the complex geological history and dramatic variations in topography [1,2,3,4]. Among the 80 known fish species documented in the Yuanjiang River, numerous species are endemic and even endangered, such as Bagarius rutilus and Poropuntius krempfi [1,2]. Furthermore, more and more invasive species such as Oreochromis niloticus, Cirrhinus mrigala and Pterygoplichthys pardalis were documented in this river during our field surveys between 2023 and 2024. Meanwhile, numerous new species have been discovered and described during last several years [5,6], indicating unexplored fish diversity are waiting to be discovered. Above all, uncovering and summarizing of the current fish diversity in the upper Red River seems to be particularly important.

Nevertheless, systematic researches of fish diversity aiming at the Yuanjiang River are rare, scattered and outdated. For example, Zhou et al. (1999) compared β diversity of three branches within the Yuanjiang River based on field surveys targeting individual branches during 1990s [7]. Liu et al. (2005) evaluated length–weight and length–length relations of twenty-three freshwater fish species in the Yuanjiang River [8]. Other relevant reports were only seen in ichthyography book and research that focused on species description [1,2]. Considering the timeliness and variability of the data on fish diversity in the Yuan River, it is highly necessary to carry out investigations and research on the fish diversity in this river.

Understanding detailed fish diversity in river is crucial for utilization and conservation [9,10,11]. To date, most taxonomic descriptions of species were largely via morphological characterizations. Nevertheless, misidentifications often occurred due to phenotypic plasticity, genotypic variation, cryptic speciation, or differing life stages [12,13,14,15,16,17,18]. DNA barcoding is a popular molecular delimitation approach for the rapid and accurate discrimination of animal species, which has demonstrated a high success rate in fish identifications regardless of the influence of life stages and/or organismal completeness [19,20,21,22,23,24]. Additionally, this technique is widely used for discovering cryptic species and novel species [25,26,27,28,29,30]. As a consequence, DNA barcoding not only enables accurate species delimitation, but also detects morphologically cryptic species and/or aids for searching new species.

Based on four field surveys during 2023 and 2024, we used DNA barcoding technique to reveal fish diversity in the Yuanjiang River. Our main goals were (1) to establish a DNA barcoding library for the ichthyofauna in the Yuanjiang River, (2) to shed new insights into fish diversity via sequence-based delimitation methods, and (3) to reveal uncertainty in species observed between DNA barcodes and morphological features.

2. Methods and Materials

2.1. Sample Collection

A total of 764 fish specimens including 64 morphological species were collected from the Yuanjiang River between 2023 and 2024 (Figure 1 and Table S1). The sampling map was generated using the ArcGIS 10.2 and modified in Microsoft Office 2017. Specimens were delimited to species following a series of monographs. Fin clips or muscle tissues were clipped and conserved in absolute ethyl alcohol for DNA extraction in the laboratory. All of the tissue samples used for genomic DNA extraction and samples were both deposited in the Freshwater Fish Museum at the Pearl River Fisheries Research Institute, Chinese Academy of Fishery Sciences.

2.2. Genomic Extraction and Sequencing

Total genomic DNA was extracted using Axygen DNA Extaction Kit following its own procedure. We then checked DNA quality using 1 % agarose gel electrophoresis. The standard barcoding fragment, ~648 bp from the 5' end of the mitochondrial COI gene, was selected for DNA barcoding marker. We utilized the fish universal primers FishF1/FishF2 and FishR1/FishR2 to amplify and sequence the target loci [31].

The polymerase chain reaction (PCR) contained approximately 100 ng of template DNA, 1μl of each primer(10 pmol), 3μl of 10× reaction buffer, 1.5μl of dNTPs (2.5 mM each), and 2.0 U of Taq DNA polymerase in a total volume of 30μl. The PCR conditions for the two fish universal primers contained initial denaturation at 95 ℃ for 5 min, 30 cycles at 94℃ for 30 seconds, 54 ℃ for 30 seconds, and 72 ℃ for 1 min, and a final extension at 72℃ for 10 min. The PCR products with high quality were sequenced bi-directionally on an ABI 3730XL DNA system (Perkin-Elmer Applied Biosystems, Foster City, USA) to decrease the occurrence of sequencing mistakes. We checked and assembled the contigs using the SEQMEN in Lasergene package (DNASTAR, Inc., WI, USA). Sequences were aligned and trimmed using the MEGA version X [32]. The aligned sequences were submitted to GenBank database with accession numbers of PV878687–PV879450 (Table S1).

2.3. Genetic Distance Calculations

We firstly calculated Kimura 2-parameter (K2P) [33] pairwise genetic distances within and among species. Secondly, to conduct a barcoding gap analysis, species-level comparisons between maximum intraspecific genetic distances and minimum distances to the nearest neighbour were calculated [34]. Lastly, we built a neighbour-joining (NJ) tree with the “pairwise deletion” option using K2P genetic distances and one thousand bootstrap replicates. The abovementioned analyses were implemented in MEGA version X.

2.4. Sequence-Based Species Delimitation

We used three species delimitation approaches to group specimens into molecular operational taxonomic units (MOTUs). Firstly, Automatic barcode gap discovery analysis was performed on the web interface (https://bioinfo.mnhn.fr/abi/public/abgd/abgdold.html) using the default value for the relative gap width (X = 1.5) and K2P genetic distance [35]. Secondly, Poisson tree processes (PTP) model in its multiple rates version was done in the web server (https://mcmc-mptp.h-its.org/; [36,37]) with a maximum likelihood (ML) tree constructed in the software RAxML-VI-HPC [38]. GTR + I + R model selected in MrModelTEST [39] was employed for ML analyses. Thirdly, Generalized Mixed Yule Coalescent (GMYC) analyses [40] using single- and multiple-threshold were also conducted using a web server (http://species.h-its.org/gmyc/). A fully resolved ultrametric gene tree as input for the GMYC analysis was reconstructed a Bayesian tree using BEAST version 1.8.2 [41]. A Yule pure-birth model and GTR + I + R substitution model were used in Bayesian reconstructions with haplotype sequences. We chose an uncorrelated relaxed lognormal clock model with the rate estimated from the data and ucld-mean parameters with a uniform prior to value of 0 as a lower boundary and 10 as an upper boundary. The remaining settings were set as the default parameters. Fifty million generations with 1000 sampling frequencies were used in the Monte Carlo Markov Chain (MCMC) run. Effective sample sizes (ESS) for all parameters were evaluated using Tracer v1.4 [42]. The final trees were summarized in a Maximum credibility tree using TreeAnnotator version 1.8.2 via discarding 50% of the total trees [41]. A majority-rule consensus out of the abovementioned three algorithms were adopted for delimitation scheme.

2.5. Cryptic Diversity

With respect to morphological species owned multiple clades and/or MOTUs with higher genetic distance, a NJ tree was built for visual inspection using K2P distances, with 1000 bootstrap replicates. In addition, mean intraspecific K2P genetic distance and genetic divergence between observed clades were also calculated. The abovementioned analyses were done using MEGA X.

3. Results

3.1. Sequence Information

We successfully obtained 764 DNA barcodes that represent 64 morphospecies, 50 genera, 15 families and 5 orders (Table S1). The number of sequences obtained per morphospecies averaged 11.94 and ranged from 1 to 42 (11 morphospecies were represented by singletons). The length barcoding sequences after aligning and trimming reached 632 base pair (bp). No deletions, insertions or stop codons were observed in any of the amplified sequences, suggesting that all of the sequences constitute functional mitochondrial COI sequences.

3.2. Genetic Distance Calculation

A hierarchical increase in the mean genetic distance was obtained from within species (0.47% ± 0.03%) to within congeners (4.12% ± 0.29%) in the K2P distance (Table 1). With respect to 53 species with more than two specimens, 49 species had an intraspecific distance less than the threshold values (2%). Nevertheless, four species, i.e., Hemiculter leucisculus, Misgurnus anguillicaudatus, Paramisgurnus dabryanus and Rhinogobius honghensis, examined relatively high intraspecific divergence (>2%). Among 53 species with more than two specimens, barcode gap analysis showed that 52 species had barcoding gaps with the minimum interspecific distance to the nearest neighbour larger than the maximum intraspecific distance (Figure 2). A regression analysis revealed that the intraspecific sequence divergence is not significantly correlated with the sample size (Spearman correlation analysis: R = 0.177, P = 0.206).

The NJ trees (Figure 3) yielded 68 clades from 64 morphospecies (including the singleton species), suggesting that several morphospecies generated multiple clusters. Four species, i.e., Hemiculter leucisculus, Opsariichthys bidens, Misgurnus anguillicaudatus and Paramisgurnus dabryanus, yielded two clusters from the NJ topology, respectively. In addition, a visual inspection of the NJ tree indicated that 60 analyzed species displayed concordant clusters of DNA barcodes.

3.3. Sequence-Based Species Delimitation

Three delimitation approaches yielded discordant results (Figure 4). The ABGD analyses conducted with the K2P models resolved 65 MOTUs. Among the 65 MOTUs, we detected that Misgurnus anguillicaudatus and Opsariichthys bidens generated two MOTUs, respectively. In addition, the ABGD analyses supported two Schistura species named S. calichromus and S. caudofurca shared a MOTUs (Figure 4). The remaining MOTUs were found to be concordant with morphological identification (Figure 4).

mPTP obtained 52 MOTUs and the number of MOTU was less than that yielded by ABGD analyses (Figure 4). Only 45 MOTUs were consistent with morphological results. mPTP analyses yielded conservative outcomes than other approaches. For example, mPTP clustered 5 morphologically distinct Schistura species into one MOTUs. In addition, six species belonging to five genera, i.e. Pseudohemiculter dispar, Anabarilius sp., Culter dabryi, Culter Alburnus, Toxabramis houdemeri and Hemicculter leuciclus into one MOTUs (Figure 4). Similar pattern was also seen in three species belonging to three genera, i.e., Discogobio yunnanensis, Garra imberba and Placocheilus caudofasciatus (Figure 4).

With respect to GMYC analyses, a total of 82 MOTUs were observed. GMYC analyses showed 11 morphological species owned 2–4 MOTUs (Figure 4), which triggered more MOTUs than the other two approaches.

The final consensus, established by a majority rule consensus, obtained 65 MOTUs. In sum, a total of 60 species were unambiguously identified by the delimitation analyses, two species (Opsariichthys bidens and Misgurnus anguillicaudatus) displayed two MOTUs and two closely related species (Schistura calichromus and S. caudofurca) were found to be displayed mixed genealogies with one MOTUs (Figure 4).

3.4. Cryptic Diversity

Two clades were detected in four species, i.e., Hemiculter leucisculus, Opsariichthys bidens, Paramisgurnus dabryanus, and Misgurnus anguillicaudatus. We found that the mean intraspecific distance ranged from 1.28% (Opsariichthys bidens) to 4.08% (Misgurnus anguillicaudatus). Furthermore, higher genetic divergence between the two clades of each species was examined, with values ranging from 2.70% (Paramisgurnus dabryanus) to 10.08% (Misgurnus anguillicaudatus) (Figure 5).

4. Discussion

The present study provides the first comprehensive DNA barcode reference library for fish species in the Yuanjiang River section in China, which includes 64 morphospecies (ca. 80% of known reported species in this river [2]) and covers most river sections including mainstem and tributaries (Figure 1). In our study, a total of 60 morphospecies (93.75% of the analyzed species) can be unambiguously delimited by DNA barcoding, providing a straightforward identification system for fish species in this river. In addition, the high intraspecific divergences observed for two species and failure delimitation detected for two closely related species demonstrate the need for further taxonomic research on these species. In sum, our findings provide robust evidence that DNA barcoding serves as an effective tool for refining extant taxonomic classifications and uncovering cryptic diversity within the ichthyofauna of the Yuanjiang River.

4.1. DNA Barcode Reference Library of Fish Species in the Yuanjiang River

The corrected mean intraspecific distance of 0.47% (SE= 0.0003) calculated for the Yuanjiang River fishes is slightly greater than the values reported in the several previous fish barcoding studies involving rivers in China, such as the Nujiang River of China section (0.41%), the middle of the Yangtze River (0.36%) and the middle and lower Pearl River (0.32%) and the Irtysh River (0.30%) [43,44,45,46]. The relatively higher level of intraspecific divergence observed in the Yuanjiang River can be explained by several species have high intraspecific genetic distance among several species. We found that the mean intraspecific genetic distances of 15 morphospecies (28.30% of analyzed species with two specimens) exceed 1.00%.

Barcoding gap analysis including 53 species with more than two samples revealed that the mean nearest neighbor distances (12.38%) were ~26-fold higher on average than intraspecific distances (0.47%), indicating that the present library is useful for the identification of the Yuanjiang River fishes. Furthermore, we demonstrated that a barcode gap was present in 52 species, with the exceptions of Schistura calichromus (Figure 2). In this study, we discovered that the interspecific divergence (3.77%) between S. calichromus and S. caudofurca was slightly lower than the maximum intraspecific divergence (3.94%). Furthermore, though NJ tree obtained monophyletic clades of S. calichromus and S. caudofurca, Bayesian trees showed the two species were mixed into one phylogenetic clade (Figure 3; Figure 4). These findings were indicative of that the two species were recent speciation and had close phylogenetic relationship. In order to disentangle the relationships of the two species, detailed studies with population-level analyses with larger sample sizes should be carried out in the future.

4.2. Cryptic Diversity Pattern

The most salient finding was that four species (Hemiculter leucisculus, Opsariichthys bidens, Paramisgurnus dabryanus and Misgurnus anguillicaudatus) generated two clades from the NJ topology, which suggested undetected genetic diversity occurring in these species. Hemiculter leucisculus is a widespread species in East Asia and is demonstrated to own undetected diversity in different drainages. For instance, Chen et al. (2017; 2021; 2022) found that H. leucisculus populations in the Pearl River generated two deep mitochondrial clades with higher genetic distance [18,47,48]. The upper Yuanjiang River is close with two independent rivers, i.e., the Pearl River and the Lancang River. Furthermore, the lower Yuanjiang River section in the Vietnam is nearby the Hainan Island drainages and they can be connected due to sea level dropped triggered by glaciation during the Late Pleistocene [49,50,51]. Collectively, two H. leucisculus clades detected in the Yuanjiang River was considered a normal phenomenon. Similar situation also can be used to explained the pattern in Opsariichthys bidens populations. In addition, Opsariichthys bidens is small and/or benthic taxa that was corroborated to harbor an unexpectedly high degree of genetic diversity within the same drainage due to limited dispersal potential, geographic isolation and geological events [52,53,54].

For Paramisgurnus dabryanus, the pattern of forming two lineages could be explained by their introductions from different sources. This species is an important economic taxon raised in many regions, including southwestern China. Introduction from various genetic resources can trigger multiple lineages. Moreover, hybridization with other species might lead P. dabryanus to generating multiple lineages due to the introduction of the mitochondrial genome from other species [55]. Similar cases can be seen in other rivers in China, such as the Nujiang River and the lower Pearl River [43,45]. Two deep clades observed in the Pond loach, Misgurnus anguillicaudatus, was an expected phenomenon for an independent river. For example, Zhong et al. (2019) argued that Misgurnus anguillicaudatus populations both in the Pearl River and the Yangtze River had four mitochondrial clades with high genetic divergence among clades [56].

4.3. Implications for the Yuanjiang River Ichthyofauna and Future Prospects

Our study provides the comprehensive geographical, taxonomic and molecular sampling of the Yuanjiang River to date, with several findings contradicting current systematic and taxonomic knowledge. The most challenging cases were those involving shallow interspecific divergence between two Schistura species and the detection of multiple, and highly divergent lineages in four species. Given that closely relationships between the two Schistura species may be recent divergence, additional molecular markers with higher substitution rates such as control region or nuclear genome dataset, would certainly help in better distinguishing cases of recent divergence. Furthermore, multiple lineages/MOTUs observed in four species can be explained by overlooked, subtle morphological differences, and cryptic diversity and unrecognized speciation events might also be responsible for this result [47]. However, the taxonomic status of these cryptic lineages delimited by mitochondrial barcoding marker is unrobust and awaits further assessments of their morphological divergence and phylogenetic placement using additional nuclear loci.

In this study, five morphologically different species belonging to four genera named as Anabarilius, Schistura, Rhinogobius and Mastacembelus cannot be assigned into exact species on the basis of current documents. This observation implies that many species in the Yuanjiang River are not described due to rare and scattered researches on this river. In addition, the Yuanjiang River ichthyofauna combinates ichthyofauna from Vietnam and Southwest China, suggesting complex fish composition. Collectively, more field surveys and taxonomic researches involving in the Yuanjiang River ichthyofauna should be conducted.

5. Conclusions

In the current study, based on high density sampling, we build a DNA barcode reference library including 64 fish species in the Yuanjiang River, China. A relatively high success rate of identification (93.75%) demonstrates the effectiveness of DNA barcoding for automated identification of the fishes in the Yuanjiang River. Newly detected mitochondrial clades in four species calls for a substantial effort to better characterize the diversity of the groups. Nevertheless, we also confirmed the limitation of DNA barcoding for automated identification of closely related species. Furthermore, five taxa cannot be delimited into species level based on morphological characteristics, indicating more taxonomic studies should be conducted for fishes in this river. Given that ca. 20% species were lack in the current study, unsampled species are waiting for updating in order to construct more comprehensive our barcode library.