DNA Barcoding of Fish Species Diversity in Guizhou, China

Abstract

Guizhou is an important ecological barrier in the upper reaches of the Yangtze River and the Pearl River basins with abundant fish species. However, fish from these regions are threatened by anthropogenic activities, including overfishing and habitat destruction. Here, we assessed the fish diversity including more than half of the species from the region using DNA barcoding (partial sequence of cytochrome c oxidase subunit I (COI) gene). We obtained 800 mitochondrial COI barcode sequences from 82 genera, 18 families and 8 orders of fishes. The average Kimura two-parameter (K2P) distances within species and genera were 0.35% and 5.44%, respectively. The average interspecific distance was 15.54 times higher than the mean intraspecific distance. Moreover, DNA barcodes revealed 175 operational taxonomic units (OTUs) based on consensus demarcation schemes. Barcoding gaps were detected in 94.81% of morphospecies. Three fish species (Schistura fasciolata, Vanmanenia pingchowensis, and Misgurnus dabryanus) have considerable intraspecific variability, and each was divided into multiple molecular operational taxonomic units (MOTUs) using molecular definition methods (Automatic Barcode Gap Discovery, Refined Single Linkage, General Mixed Yule Coalescent, and Poisson Tree Processes), possibly indicating the occurrence of cryptic species. Altogether, our study reveals the complex diversity of fish species in Guizhou Province, serving as a reference for the conservation and monitoring of fish species in this region.

1. Introduction

Fish are an important part of biodiversity in the composition of animal groups [1]. Because of their wide distribution and rich species diversity, fish not only have great nutritional value but also play an important role in ecological conservation [2]. However, human activities such as overfishing, sewage discharge, and river dam building have adversely affected the fish’s survival recently [3]. Therefore, the classification and identification of fish species are basic requirements for diversity conservation.

Past studies on the identification and classification of fish species relied on morphometric and taxonomic characteristics, such as the number of fin rays and lateral scales, body size, and color [1]. However, the morphological characteristics of fish species are not always consistent in different stages of development [4,5]. Morphologically incomplete and hidden species of fish cannot be identified accurately, which will make the identification more difficult and even result in incorrect identification. Currently, the lack of experienced fish taxonomists has led to a knowledge gap in the species-level identification of some groups. Moreover, there are also differences in the classification description of the same fish in different literature, which can lead to misidentification by novices. Thus, molecular techniques have been widely applied to address the ambiguities.

DNA barcoding has now been developed to serve as an efficient and fast identification of species [6]. DNA barcode technology can identify animal species using a partial sequence of the mitochondrial cytochrome c oxidase subunit I (COI) [7]. Moreover, DNA barcoding identification of species with extremely similar morphology, species at different stages of growth, and species that are morphologically incomplete is more accurate than morphological identification [8,9,10]. Since the DNA barcoding technique was applied by Ward et al. [11] to identify fish species successfully, fish DNA barcoding has been adopted globally, and the technology has been widely used [12,13,14,15,16]. The COI gene is the genetic marker most utilized in fish DNA barcoding research because it can accurately identify 93% of freshwater fish [17]. Nevertheless, the COI gene has some disadvantages of mitochondrial genes, such as incomplete lineage sorting phenomena and gene introgression, which may lead to misidentification [18,19]. Thus, the barcoding method is not aimed at replacing taxonomy but rather at providing a useful additional tool in case of unclear identification (in particular for disputed fish species) to assess biodiversity. As a result, various analysis methods combined with traditional morphological identification are used to define species in DNA barcoding studies to improve the accuracy of taxonomic identification, enhance the discovery of cryptic species and enrich the genetic diversity of species [20,21,22,23,24,25].

Guizhou is located at the confluence of the upper reaches of the Yangtze River and Pearl River basins. It is also the center of a karst region in southwest China, with developed karst landforms and numerous underground rivers [26,27]. The complex topography and hydrology provide the region with a wide variety of habitats and rich biodiversity [28,29]. Three comprehensive and systematic fish resource surveys have been conducted based on morphology in Guizhou Province in the past 30 years [30,31,32]. The latest checklist of freshwater fishes of Guizhou contains 288 species [32]. However, human activities in the freshwater environment of Guizhou, such as the excessive construction of water conservancy facilities, sewage discharge and overfishing, have escalated in recent years, leading to the loss of freshwater habitats and serious damage to the structure of fish resources [33]. As a result, a large-scale taxonomic assessment of ichthyofauna needed to be conducted for further conservation and sustainable management of fish resources in the region.

The present study reappraises Guizhou ichthyofauna diversity through a DNA barcode–based approach to species delimitation. This study substantially improves our understanding of genetic diversity of the river systems of the region. Additionally, the DNA barcodes generated in this study will be available for monitoring and conservation of fish diversity.

2. Materials and Methods

2.1. Specimen Collection and Morphological Identification

Field sampling was conducted in rivers, lakes, and reservoirs in Guizhou, China, from May 2017 to January 2022. The sampling sites (N = 103) of two river basins in Guizhou are shown in Figure 1. Fishing gear, such as traps, fishing rods, and gillnets, were used to catch the fish. Some fish species were obtained from fishermen. The specimens were preserved in 75% ethanol, then transported to the laboratory. Voucher specimens were deposited in the School of Life Sciences, Guizhou Normal University.

All specimens were identified based on morphological characteristics as described in Wu [30] and Yang et al. [32]. Valid names of the fish were checked using Fishbase (https://www.fishbase.de/, accessed on 28 May 2022) and Eschmeyer’s Catalog of Fish (https://researcharchive.calacademy.org/research/ichthyology/catalog/fishcatmain.asp, accessed on 28 May 2022).

2.2. DNA Extraction, PCR Amplification, and Sequencing

Genomic DNA was extracted from the muscle tissue using a standard high-salt method as previously described [34]. A 1% agarose gel electrophoresis was then used to confirm the integrity of the genomic DNA, then stored at −20 °C for further use. The COI gene sequences of samples were amplified using the universal primers described by Ward et al. [11]: FishF1- TCAACCAACCACAAAGACATTGGCAC and FishR1- TAGACTTCTGGGTGGCCAAAGAATCA. The polymerase chain reaction (PCR) mixture contained 1μL of template DNA (50–100 ng/μL), 1 μL of each primer (10 μM), 17.5 μL of 2xTaq Plus MasterMix (CoWin Biosciences, Beijing, China) (total volume 35 μL). The PCR conditions were: initial denaturation at 95 °C for 5 min; 35 cycles of denaturation at 95 °C for 1 min, annealing at 55 °C for 30 s and extension at 72 °C for 1 min, and a final extension at 72 °C for 10 min. The following specific COI primers for Gambusia affinis were designed based on the mitochondrial genome of G. affinis (GenBank accession number: NC_004388): GA_F-TCAACCAATCATAAAGATATCGGCAC and GA_R- TAAACTTCCGGATGCCCGAAAAACCAG. PCR cycling was conducted as described for the universal primers. The following specific COI primers for Micropercops swinhonis were designed based on the mitochondrial genome of M. swinhonis (GenBank accession number: NC_021763): MS_F- TCAACTAATCATAAAGACATTGGTAC and MS_R- TATACTTCTGGGTGACCAAAAAACCAA. Except for the annealing temperature at 47 °C, the other PCR cycling conditions were similar to those for the universal primers.

The amplified DNA was fractionated via electrophoresis on 1% agarose gels. The lengths of fragments were determined by comparing them against the DL2000 DNA marker (TaKaRa, Beijing, China). The PCR products were taken to Sangon Biotech. Co., Ltd. (Shanghai, China), for purifying and bidirectional sequencing. Sequence analysis of the PCR products was carried out by Sanger sequencing using the BigDye Terminator v3.1 Cycle Sequencing kit on a 3730xl DNA analyzer (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions.

2.3. Data Analysis

The sequences were manually assembled using Seqman in the DNASTAR Lasergene package (DNASTAR, Inc. Madison, WI, USA). MEGA version 7.0 [35] was used to align and edit all sequences using MUSCLE (Multiple Sequence Comparison by Log-Expectation) default parameters. Identification Systems (IDs) and the Basic Local Analysis Search Tool (BLAST) from Barcode of Life Data Systems (BOLD) and the National Center for Biotechnology Information (NCBI), respectively, were used to compare the sequences to the standard sequence. Reassessment was conducted for some samples that were difficult to identify and those that were not successfully compared in the database. The assessment criteria were based on the distribution location of the species, the species morphology, and the position of the sequences on the phylogenetic tree [23]. The R package ape 5.0 [36] was used to calculate Kimura two-parameter (K2P) [37] pairwise genetic distances. The pairwise K2P genetic distances matrix and the R package Spider 1.5 [38] were used to determine the maximum intraspecific and nearest-neighbor genetic distances. The phylogenetic tree was reconstructed based on the Bayesian inference (BI) and Maximum Likelihood (ML) methodology using MrBayes 3.2.7 [39] and RAxML [40], respectively. The GTR+I+G model was selected as the best model using MrModeltest 2.4 [41]. The COI sequences were subsequently uploaded to the BOLD online system as part of a project entitled “Freshwater Fish in Guizhou Province”.

Four methods were used to determine molecular operational taxonomic units (MOTUs) from the DNA barcodes for the estimation of genetic diversity using single mitochondrial gene sequences: (a) Refined Single Linkage (RESL) analysis was used to generate OTUs from sequences based on the “Cluster Sequences” in the BOLD [42]; (b) Automatic Barcode Gap Discovery (ABGD) analysis was conducted using the web interface (https://bioinfo.mnhn.fr/abi/public/abgd/, accessed on 28 May 2022), Kimura (K80) TS/TV distance with X (relative gap width) of 1.5 and Nb bins (for distance distribution) of 20 [43]; (c) Bayesian implementation of the Poisson Tree Processes (bPTP) was used to reconstruct a maximum likelihood (ML) tree using the RAxML server with set 1000 bootstrap replicates. The obtained trees were visualized using FIGTREE v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/) and then it was conducted on the web server (http://species.h-its.org/) [44]. (d) General Mixed Yule Coalescent (GMYC) analysis was first run using the Bayesian Ultrametric tree in BEAST.v2.6.4 (http://www.beast2.org/, accessed on 28 May 2022) with the following settings: Yule model, strict clock, GTR+I+G model, Monte Carlo Markov chains (MCMC) for 50 million generations [45]. The Effective sample size (ESS) values were evaluated using TRACER.v1.7.2 (http://beast.community/tracer, accessed on 28 May 2022) with the previously used parameters. A maximum clade credibility tree was then constructed using TreeAnnotator.v2.6.4. Finally, the tree was visualized using R packages ape and splits [46]. A final delimitation scheme was established based on a majority-rule consensus among the four delimitation analyses and one morphological identification performed.

3. Results

A total of 159 morphospecies belonging to 82 genera, 18 families and 8 orders were identified (Table S1). Twenty-four species were not identified at the species level. Moreover, 93 and 66 species were collected in the Pearl River basin and Yangtze River basin, respectively. The mean number of specimens per species was 9 (1–30). Twenty-five species had only a single sequence.

A total of 800 barcode sequences were successfully obtained from the samples. All of the sequences formed 651 bp after editing, and no stop codons, insertions, or deletions were found in any of the sequences. The intraspecific genetic distances ranged from 0 to 3.14% (average; 0.35%). Macropodus opercularis, Misgurnus dabryanus, Vanmanenia pingchowensis, and Schistura fasciolata had a genetic distance of more than 2% (Table S3). The interspecific genetic distance ranged from 0.65% to 11.02% (average; 5.44%) (15.54 times that of intraspecific genetic distances) (

Table 1. Summary of K2P distances within species and genus levels of fishes.

Level	Min Distance (%)	Mean Distance (%)	Max Distance (%)	SE (%)
Species	0	0.35	3.14	0.10
Genus	0.65	5.44	11.02	0.64

). The pairwise K2P genetic distances analysis revealed that five species (Discogobio brachyphysallidos, Hemibarbus maculatus, Hemibarbus medius, Paraqianlabeo lineatus, Sinocyclocheilus cyphotergous, and Sinocyclocheilus multipunctatus) did not have a barcode gap (Figure 2, Table S3).

The MOTU delimitation analyses identified different numbers of MOTUs (Figure 3, Table S2). For instance, RESL, ABGD, GMYC, and bPTP identified 166, 161, 185, and 200 delimited MOTUs, respectively. The number of results based on molecular delimitation was higher than that of morphological identification, especially with the bPTP method. The final consensus obtained 175 OTUs. The distribution of maximum intraspecific distance and the nearest-neighbor genetic distances for OTUs did not overlap and showed a significant barcoding gap (Figure 2). Thirteen morphospecies (M. swinhonis, Mastacembelus armatus, M. opercularis, Pterocryptis anomala, V. pingchowensis, S. fasciolata, Microphysogobio elongatus, Acrossocheilus yunnanensis, D. brachyphysallidos, Parasinilabeo assimilis, Pseudocrossocheilus bamaensis, H. maculatus, M. dabryanus) were split into multi-OTUs based on final consensus, while Tachysurus brachyrhabdion and Tachysurus gracilis, S. multipunctatus and S. cyphotergous, Acrossocheilus longipinnis and Acrossocheilus iridescens, H. medius and Hemibarbus labeo were recovered under single OTUs (Table S2).

Moreover, S. fasciolata, V. pingchowensis, and M. dabryanus had high intraspecies variation and more than one MOTU based on the four delimitation methods (Table S2), indicating that they might be cryptic species.

Seventeen specimens of the genus Schistura were also collected, of which eight were identified as S. fasciolata. However, nine were identified at the species level and were distributed into four morphologically distinct clusters (Schistura sp.1, Schistura sp.2, Schistura sp.3, and Schistura sp.4). The four clusters had an average interspecific genetic distance of 3.80%. Additionally, the four clusters were divided into a single MOTU using RESL, GMYC, and bPTP. The molecular delimitation results were similar to morphological identification results. However, the average intraspecific genetic distance of morphologically identified S. fasciolata was 2.7% (0 to 4.32%), above the estimated delimitation threshold (2%). The phylogenetic tree shows two large clades of the genus Schistura (Figure 4). S. fasciolata had eight analyzed specimens and showed a putative cryptic species represented by five MOTUs in all delimitation methods, supported by five clades in the phylogenetic tree. Only the F170 was collected from the Yangtze River system, while the other species were collected from the Pearl River system. Schistura sp. formed four distinct clades in the tree, of which four branches were collected from different locations in the Pearl River system.

Sixteen specimens of the genus Vanmanenia were collected, of which 13 and 3 were identified as V. pingchowensis and Vanmanenia sp., respectively, based on morphology. The other three specimens were identified at the genus level (i.e., Vanmanenia sp.). The V. pingchowensis and Vanmanenia sp. have a high genetic distance (13.56%) and reached intergeneric divergence. However, V. pingchowensis also has high intraspecific variation (3.14%), divided into three OTUs using molecular methods. Sixteen specimens of the genus Vanmanenia formed two larger clades in the tree (Figure 5). Moreover, 13 specimens of V. pingchowensis from the Pearl River basin and the Yangtze River basin have two distinct clades, consistent with MOTUs results.

Eight specimens of M. dabryanus were collected from seven locations. The specimens were split into two MOTUs via delimitation methods and have high intraspecific variability (2.27%). The phylogenetic tree showed that M. dabryanus had formed two clades, not related to their geographical origin (Figure 6).

4. Discussion

To our knowledge, our work represents the first molecular survey of ichthyofauna in Guizhou, China. Specifically, 800 DNA barcode sequences were detected for 159 species, representing 55.21% of the number of fish species reported in the region [32].

Genetic distances are analyzed mainly using the barcode gap and the genetic threshold. A species has a barcode gap when the minimum interspecies distance is greater than the maximum intraspecies distance [47]. According to Hebert et al. [7], a threshold of 10 times or more interspecific and intraspecific genetic variation is sufficient for the successful application of DNA barcoding. The rule has also been applied in other DNA barcoding studies, such as Cuban freshwater fishes [13], the ichthyofauna of the Yangtze River [48], and Uzbek fishes [1], thereby further confirming that this approach is an effective way to distinguish between fish species. In this study, the average interspecific genetic distance was 15.54 times higher than the average intraspecific genetic distance. Moreover, 94.81% of the morphospecies had significant barcode gaps. These results confirmed the effective use of DNA barcoding for fish species diversity studies in Guizhou. However, the higher DNA barcoding capacity could be because most genera in the dataset had only one representative species and a low number of homologous species. Bagley et al. studied fish DNA barcoding in Cerrado headwater streams in Central Brazil and reached the same conclusion [49]. This may be common due to unclear species boundaries or limited geographic sampling of intraspecific variation [50]. Furthermore, the mean interspecific genetic distance (5.44%; 0.65%−11.02%) was lower for fish species in this study than for Hainan fish (14.7%) [51], Gansu fish (9.99%) [52] and Henan fish (14.54%) [53]. This finding could be attributed to the lower interspecific genetic variation vis-à-vis the intraspecific genetic variation in the genera studied (Hemibagrus, Hemibarbus, and Hemiculter). In addition, it also may be due to incomplete lineage sorting and recently diverged species since their molecular diversity is not robust [54].

However, the molecular delimitation methods and morphology-based identifications had some inconsistencies. These results confirm the benefits of combining several species delimitation methods. We found that multiple new OTUs of Rhinogobius, Schistura and Vanmanenia were identified only to the genus level and are now awaiting description. Furthermore, 14 morphospecies had a shallow genetic distance (<2%) and consisted of more than one MOTU based on at least two molecular delimitation methods (Table S2). This result highlights that the diversity of the freshwater fish of this region is probably largely underestimated. However, a study assessing a wider geographical area is needed to assess the cryptic diversity of the morphospecies.

The number of MOTUs split by ABGD and RESL was similar to that of morphologically identified species. However, the number of MOTUs produced by GMYC and bPTP was higher than that of morphospecies. Moreover, the bPTP approach defined the highest number of MOTUs (200 MOTUs), consistent with the findings of freshwater fishes of Indo-Myanmar biodiversity hotspot [3], Korean curved-horn moths [55] and the Congolese and Lower Guinean ichthyological provinces [56]. Tree-based species delimitation of the bPTP method successfully separated Arossocheilus iridescens, A. longipinnis, T. brachyrhabdion, and T. gracilis in a single species, consistent with morphological identification (Table S2). The 800 fish samples were divided into 175 OTUs based on the combined five species definition methods. All OTUs had significant barcode gaps, indicating that the consistent definition results were acceptable.

The interspecies had close relatives and shared a single MOTU based on multiple species (H. medius and H. labeo, A. iridescens and A. longipinnis, T. brachyrhabdion and T. gracilis, S. cyphotergous and S. multipunctatus). Although S. cyphotergous and S. multipunctatus have significantly different morphology, the two species are closely related in the evolutionary tree of Sinocyclocheilus [57]. To date, the molecular delimitation methods employed have failed to recognize fish species because of the absence of interspecific diversity in the same river basin (e.g., genera Triplophysa and Tachysurus of the Yangtze River [48], genera Schizothorax of the Salween River [58], and Carassius of the Mediterranean basin) [1]. The close relationship may be because the divergence of the species occurred recently; thus, the genetic distance has not reached a threshold to allow for separation via the molecular delimitation method. However, these species do not share any haplotype, possibly due to the absence of interspecific hybridization and gene introgression.

The phenotypic plasticity significantly influenced morphological identification, leading to high intraspecific divergence in morphospecies [3]. In this study, the four species (M. dabryanus (2.27%), S. fasciolata (2.70%), V. pingchowensis (3.14%), and M. opercularis (2.20%)) have high intraspecific variability. Each was divided into multiple MOTUs using the molecular definition method. Other studies assessing a wide geographical area have identified cryptic species indicating allopatric speciation [59]. The three cryptic species identified here also support this view. However, the level of genetic differentiation among M. dabryanus populations did not show significant geographical differentiation, inconsistent with previous reports [60]. Furthermore, two specimens of M. opercularis were collected from the same place in the Pearl River basin. Therefore, sympatric coexistence of these two specimens and fewer sample numbers may lead to a high level of intraspecific genetic diversity. In another study, M. opercularis has a low-density population, genetic diversity, and geographical differences [61]. As a result, M. opercularis was not considered a cryptic species. Future studies should collect more samples to provide an in-depth analysis of the genetic diversity of Macropodus.

5. Conclusions

In many studies, DNA barcoding proved to be an effective supplementary tool for accurate species identification and biodiversity research. In this study, a large-scale molecular evaluation of wild fish species was conducted in Guizhou Province using DNA barcoding. The results confirmed that DNA barcoding can effectively identify fish in Guizhou Province. About 94.81% of morphological species had significant genetic differences. Moreover, the number of MOTUs divided via the four molecular definition methods exceeded the number of species identified based on morphology. High intraspecific variation indicates the existence of cryptic species or hidden diversity in the fish fauna of Guizhou. In addition, the newly found unknown fishes indicates the possible existence of new species. However, future studies should increase fish samples for a broader assessment of genetic diversity, including the use of nuclear genes and other molecular markers to further verify the boundaries of these species. In summary, combining molecular delimitation methods and morphology methods can correctly identify the fish species in the region.