Large-scale DNA barcoding of the subfamily Culterinae (Cypriniformes: Xenocyprididae) in East Asia unveils a geographical scale effect, taxonomic warnings and cryptic diversity.

Geographical scale might be expected to impact significantly the efficiency of DNA barcoding as spatially comprehensive sampling provides opportunities to uncover intricate relationships among closely related species and to detect cryptic diversity for widespread taxa. Here, we present a DNA barcoding study on a Xencyprididae subfamily (Culterinae) involving the production of 998 newly generated DNA barcodes from East Asian drainages (80 localities). Together with 513 barcodes mined from BOLD and GenBank, a reference library consisting of 1511 DNA barcodes (116 localities) for 42 species was assembled, accounting for 66% of known Culterinae species. Intraspecific genetic distances are positively correlated to geographical scale, while a negative correlation is detected between interspecific genetic distances and geographical scale. The present study demonstrates that geographical scale influences the efficiency of DNA barcoding by narrowing the width of the barcoding gap. DNA-based species delimitation analyses delimited 44 molecular operational taxonomic units (MOTUs). Rampant cryptic diversity is detected within eight species with multiple MOTUs, whereas 25 species present mismatch between morphological and molecular delimitations. A total of 18 species are lumped into nine MOTUs due to low interspecific divergence and/or mixed lineages. Several MOTU divergences are hypothesized to relate to known biogeographical barriers and geological events during the Pliocene and Pleistocene. This study provides new insights into the taxonomy and phylogeography of the subfamily Culterinae.

INTRODUCTION

DNA barcoding, proposed as a solution to overcome the rarefaction of taxonomists worldwide (i.e., taxonomic impediment), provides a reliable framework for sustainable species identification and assessment of neglected or poorly investigated taxonomic groups (April et al., 2011; Chen et al., 2015; Galimberti et al., 2012; Galimberti et al., 2021; Hebert et al., 2003; Shen et al., 2019). Since this approach, and associated data quality standards, was developed, a large number of empirical DNA barcoding studies have been conducted in various taxa (April et al., 2011; Chen et al., 2015; Delrieu‐Trottin et al., 2019; Galimberti et al., 2012; Galimberti et al., 2021; Hebert et al., 2003; Hebert et al., 2004; Hubert et al., 2008; Hubert et al., 2012; Pereira et al., 2013; Shen et al., 2019; Sholihah et al., 2020). Most of these studies yielded a broad match between morphological identifications and DNA barcode clusters (Ball et al., 2005; Breman et al., 2016; Delrieu‐Trottin et al., 2020; Hebert et al., 2003; Hebert et al., 2004; Kim et al., 2020; Lara et al., 2010; McCusker et al., 2013; Panprommin et al., 2019; Ribeiro et al., 2012; Smith et al., 2005), resulting in the detection of a barcoding gap, which is the lack of overlap between the distributions of the maximum intraspecific and minimun interspecific genetic distances (Meyer & Paulay, 2005). However, most previous reports focused on particular areas, such as a nation or a biogeographical province (Ball et al., 2005; Breman et al., 2016; Hebert et al., 2003; Hebert et al., 2004; Kim et al., 2020; Lara et al., 2010; McCusker et al., 2013; Panprommin et al., 2019; Ribeiro et al., 2012; Smith et al., 2005). If only restricted areas are involved, only a subset of global diversity of resident lineages is considered, leading to two shortcomings, which artificially inflate the width of the barcoding gap: (i) in spatially structured populations, only a fraction of intraspecific genetic diversity is sampled, and maximum intraspecific genetic distances are underestimated; and (ii) if allopatric speciation mostly contributed to species pools, sister‐species are under‐represented and minimum genetic distances between species are over‐estimated (Bergsten et al., 2012). As such, regional species assemblages yield a greater match between morphological delimitations and DNA barcode clusters than entire lineages, a trend calling for broad‐scale assessments in the case of highly diversified lineages.

The Xenocyprididae subfamily Culterinae (Ostariophysi: Cypriniformes) is a group of Asian endemic fishes (Figure 1), mostly residing in East Asian drainages (Froese & Pauly, 2020; Luo & Chen, 1998). A total of 64 recognized species, belonging to 18 genera, have been documented in this group (Dai et al., 2005). Among recognized species, many are economically important in East Asia (e.g., Chanodichthys spp., Culter spp., Hemiculter spp., Parabramis spp. and Megalobrama spp.) and are targeted by inland fisheries due to their medium to large size (Froese & Pauly, 2020; Luo & Chen, 1998). Additionally, several species such as Anabarilius spp., Pogobrama barbatula, Hainania serrata, Ancherythroculter wangi and Megalobrama pellegrini are on the list of endangered species in the “Red List of China's Vertebrates” (Jiang et al., 2016). As such, this subfamily has attracted much research interest in recent decades because of their conservation status (Cao & Wang, 2010; Feng et al., 2009; Jiang et al., 2018; Li et al., 1993; Zheng, 1989).

FIGURE 1

Selected species of Culterinae, which illustrate the diversity of the subfamily. (a) Ancherythroculter lini Luo, 1994; maximum standard length (SL) of 19.1 cm. (b) Culter recurviceps (Richardson, 1846); maximum total length (TL) of 54.6 cm. (c) Culter mongolicus (Basilewsky, 1855); maximum SL of 100.0 cm. (d) Chanodichthys erythropterus (Basilewsky, 1855); maximum TL of 102.0 cm. (e) Megalobrama pellegrini (Tchang, 1930); maximum SL of 25.1 cm. (f) Toxabramis houdemeri Pellegrin, 1932; maximum SL of 14.8 cm. (g) Metzia formosae (Oshima, 1920); maximum SL of 7.3 cm. (h) Parabramis pekinensis (Basilewsky, 1855); maximum TL of 55.0 cm. Maximum SL and TL data from fishbase (Froese & Pauly, 2020)

In the subfamily Culterinae, many closely related and morphologically similar species have been described (Froese & Pauly, 2020; Luo & Chen, 1998), and several cases of discordant results between morphological and genetic species delimitation have been observed (Luo, 2016; Wang et al., 2019; Xie et al., 2012), casting doubt on whether accurate species identification can be performed using DNA barcodes for this subfamily. For example, several morphologically different species within the genera Megalobrama or Culter were found to be nested within one lineage and/or share mitochondrial and nuclear haplotypes (Luo, 2016; Wang et al., 2019; Xie et al., 2012). Moreover, crossing experiments among the four Megalobrama species showed that they can hybridize (Zhang et al., 2014), introducing further difficulties in accurate molecular species identification.

Given some widespread Culterinae species are distributed across East Asia, a proportion of these species appear to form multiple hidden lineages under the influence of biogeographical barriers along with palaeogeological events. The most famous biogeographical barrier for freshwater fishes in China is the Nanling Mountains, which separates the Yangtze and the Pearl Rivers (Figure 2). This barrier has been demonstrated to have shaped distinct ichthyofauna in these two watersheds (Chen et al., 1986; Zheng, 1989) and to have led to a substantial number of cryptic lineages during the uplift of the Qinghai–Tibetan plateau and the evolution of the Asian monsoon (An et al., 2001; Chen et al., 2017; Perdices et al., 2005; Yang et al., 2009; Zhong et al., 2019). Other biogeographical barriers such as the Qiongzhou Strait and the Yunkai Mountains were demonstrated to have accelerated the divergence of local fish taxa in southern China (Chen et al., 2007; Chiang et al., 2013; Yang et al., 2012; Yang et al., 2016). Accordingly, if genetically structured populations are observed across those barriers, inflated intraspecific genetic divergences are expected for widespread taxa.

FIGURE 2

Sample sites for the Culterinae species included in the present study. Details of the 116 sites and collection of DNA barcoded specimens are provided in Appendix S1. Green and red circles represent sampling sites for novel obtained barcodes and downloaded barcodes, respectively. Red and yellow shadows indicate the Yangtze and the Pearl River drainage, respectively. This map was created in arcgis version 10.2 (http://www.esri.com/arcgis/about‐arcgis)

Here, we present the results of a large‐scale DNA barcoding campaign for Culterinae in East Asia. A total of 1511 DNA barcodes from 116 locations, covering most important rivers in East Asia, were included. Furthermore, widely distributed species were sampled across different biogeographical regions, providing suitable candidates to probe biogeographical patterns. Thus, the goals of the present study were: (i) to assess the performance of DNA barcoding in capturing Culterinae species boundaries in large‐scale sampling, (ii) to examine the impact of population genetic structure and geographical distance on the relative distribution of intraspecific and interspecific genetic distances, and (iii) to estimate the proportion of cryptic diversity and discuss its origin in the light of the biogeographical history of the region. Finally, this study aims to streamline the inventory of Culterinae species and promote further taxonomic studies.

MATERIALS AND METHODS

Sampling and identification

A total of 998 Culterinae specimens were collected between April 2010 and July 2020 across 80 sites that covered the Yangtze River, the Pearl River, the Yellow River, Heilongjiang River, the Huaihe River, and many independent and small coastal rivers (Figure 2; Table S1). Voucher specimens were morphologically identified using the monograph “Fauna Sinica Osteichthyes, Cypriniformes II, Ichthyography of the Pearl River and Freshwater Fishes of Guangxi, China” (Luo & Chen, 1998; Zheng, 1989; Zhou & Zhang, 2005). A fin clip or a muscle biopsy was taken from each specimen and fixed in 99% ethanol solution for further genetic analyses. Tissues and voucher specimens were deposited at the Pearl River Fisheries Research Institute, Chinese Academy of Fishery Science.

DNA sequencing, data preprocessing and assembly

Genomic DNA was extracted using the Genomic DNA Isolation Kit (Axygen) according to the manufacturer's instructions. The standard 652‐bp fragment from the 5′ end of the mitochondrial cytochrome c oxidase subunit I gene (COI) was amplified using the fish universal primers FishF1 and FishR1 (Ward et al., 2005). Amplifications included an initial denaturation step at 95°C for 5 min followed by 30 cycles of denaturation at 95°C for 1 min, annealing at 50°C for 45 s, and extension at 72°C for 45 s followed by a final extension at 72°C for 10 min. All amplified fragments were purified using 1.0% low‐melting agarose gel electrophoresis and bidirectionally sequenced with the same primer pair using an ABI PRISM 3700 (Applied Biosystems) automatic DNA sequencer. COI contigs were checked and assembled using the SEQMEN procedure in the lasergene package (DNASTAR). Sequences were aligned and trimmed to the same length using mega version 6.0 (Tamura et al., 2013) and further translated into amino acids to check for the potential occurrence of stop codons. Aligned sequences, primer pairs, trace files, taxonomic information and collection data were deposited at the Barcode of Life Data system (BOLD; Ratnasingham & Hebert, 2013) in the project CUBAR “Barcodes of Culterinae species” (dx.doi.org/10.5883/DS‐CUBAR). Alongside newly generated DNA barcodes, DNA barcode records belonging to Culterinae species were mined from GenBank and BOLD, where sample sources were available (Table S1) and sequences exceeded 600 bp.

Genetic distance analyses

Using mega 6, we calculated Kimura 2‐parameter (K2P) (Kimura, 1980) pairwise genetic distances within and among species. To check for the presence of a barcoding gap, species‐level comparisons between maximum intraspecific genetic distances and minimum distances to the nearest neighbour (Meyer & Paulay, 2005) were performed in mega 6. To explore the effect of geographical scale on genetic distance estimates, individuals from multiple geographical locations were included. Linear regression and ANOVAs were performed using spss 18.0 to assess the relationship between species sample size and maximum intraspecific divergence, and between maximum spatial distance and maximum intraspecific divergence. The effect of geographical scale on the genetic distance to the closest species was investigated by creating five geographical distance categories, <1, <50, <200, <1000 and <5000 km. For each category, all interspecific genetic distances were calculated, and the minimum was recorded for each species. Finally, a neighbour‐joining (NJ) tree was reconstructed for a visual inspection of the data, using K2P genetic distances with the “pairwise deletion” option in mega 6. We used 1000 bootstrap replicates to assess the branch support of the NJ tree.

Sequence‐based species delimitation

To group specimens into molecular operational taxonomic units (MOTUs), defined as diagnosable molecular lineages (Avise, 2012; Moritz, 1994), four sequence‐based species delimitation methods were used, including Refined Single Linkage (RESL) (Ratnasingham & Hebert, 2013), Automatic Barcode Gap Discovery (ABGD) (Puillandre et al., 2012), Poisson Tree Process (PTP) (Kapli et al., 2017; Zhang et al., 2013), and the General Mixed Yule‐coalescent (GMYC) (Fujisawa & Barraclough, 2013). A general delimitation scheme was established based a majority‐rule consensus of the four algorithms used. For the sake of clarity, species identified by morphological features are referred to as species, and species delimited using DNA barcodes are referred to as MOTUs.

RESL analysis was used to produce Barcode Index Numbers (BINs) as implemented in BOLD. ABGD analysis was performed on the web interface (www.abi.snv.jussieu.fr/public/abgd/) using the default value for the relative gap width (X = 1.5) and K2P genetic distance. PTP in its multiple rates version was run using the web server (http://mptp.h‐its.org). As the mPTP algorithm needs a phylogenetic tree as an input file, a maximum‐likelihood (ML) tree was built using raxml‐vi‐hpc (Stamatakis, 2006). ML analyses were conducted using the GTR + I + Γ model, according to the model selection procedure implemented in mrmodeltest (Nylander, 2004). Single‐ and multiple‐threshold GMYC analyses were conducted using a web server (http://species.h‐its.org/gmyc/). As GMYC requires a fully resolved, ultrametric gene tree as input for the analysis, we reconstructed a Bayesian tree using beast version 1.8.1 (Drummond & Rambaut, 2007). Bayesian reconstructions were conducted using a Yule pure‐birth model and GTR + I + Γ substitution model. Duplicated sequences were collapsed into haplotypes using alter (Glez‐Pena et al., 2010) prior to running beast analyses. An uncorrelated relaxed lognormal clock model was used with the rate estimated from the data and ucld‐mean parameters with a uniform prior value of 0 as a lower boundary and 10 as an upper boundary. All other settings were kept as the default values. Markov chain Monte Carlo (MCMC) runs were performed using 50 million generations, sampled every 1000 trees. Stability around unimodal posteriors and effective sample sizes (ESS) for all parameters were assessed using tracer version 1.5 (Rambaut & Drummond, 2007). After removing 50% of the trees as burn‐in, trees were summarized in a maximum credibility tree using treeannotator version 1.8.1 (Drummond & Rambaut, 2007).

Cryptic diversity and phylogeography

For any species pairs or groups of species with shallow genetic divergence, a network approach based on the median‐joining algorithm was used to explore the relationships among haplotypes, as implemented in network 4.6 (Bandelt et al., 1999).

With respect to species dissolving into multiple MOTUs, an NJ tree was built for visual inspection using K2P distances, with 1000 bootstrap replicates, and lineages were projected onto a map. In addition, divergence times among MOTUs were estimated using beast 1.8.1 following the procedure mentioned previously. Trees were reconstructed using a strict‐clock model, a substitution rate of 1.2% per million years (Bermingham et al., 1997), a Yule diversification model and GTR + I + Γ substitution model. Analyses were conducted with 100 million MCMC generations, sampled every 1000 generations. Trees were summarized with treeannotator version 1.8, and the first 50% of sampled trees were removed.

RESULTS

Within 998 newly acquired samples between 2010 and 2020, 26 Culterinae species, belonging to 13 genera, were delimited and identified based on their external morphology. Among them, 998 DNA barcodes were successfully produced. All the sequences were above 600 bp in length without stop codons or insertions, indicating collected sequences represent functional coding regions. Alongside newly generated DNA barcodes, 513 DNA barcodes belonging to 35 Culterinae species were mined from GenBank and BOLD (Table S1). After aligning newly generated and mined sequences, the final alignment consisted of 602 bp for 1511 specimens from 116 locations, covering most drainages in East Asia (Figure 2), and corresponding to 42 species and 15 genera (Table S1). The number of sequences per species ranged between one and 261, with an average of 34, originating from two to 38 sampling sites (Table 1).

TABLE 1

Summary statistics of the number of individuals per species (N), the number of sampling sites per species (N P), maximum intraspecific genetic distance (D M), genetic distance to the closest neighbour (D N) and the geographical extent of sampling in kilometres (DIS)

Species	N	N P	D M	D N	DIS
Anabarilius brevianalis	1	1	NA	0.0468	NA
Anabarilius grahami	3	1	0.0033	0.0698	0
Anabarilius liui chenghaiensis	2	1	0.000	0.0101	0
Anabarilius liui yalongensis	1	1	NA	0.0101	NA
Ancherythroculter kurematsui	2	1	0.0000	0.0017	0
Ancherythroculter lini	38	12	0.0050	0.0441	713
Ancherythroculter nigrocauda	15	3	0.0017	0.0017	536
Ancherythroculter wangi	1	1	NA	0.0017	NA
Chanodichthys compressocorpus	1	1	NA	0.0325	NA
Chanodichthys erythropterus	103	20	0.0707	0.0325	3526
Culter alburnus	133	28	0.0118	0.0113	1646
Culter dabryi	23	7	0.0033	0.0017	2829
Culter mongolicus	54	14	0.0396	0.0196	3560
Culter oxycephaloides	1	1	NA	0.0311	0
Culter recurviceps	261	29	0.0169	0.0113	877
Hemiculter bleekeri	70	9	0.0101	0.0036	2911
Hemiculter leucisculus	220	38	0.0556	0.0266	3657
Hemiculter lucidus	6	1	0.0050	0.0036	0
Hemiculter tchangi	20	3	0.0067	0.0332	136
Hemiculterella sauvagei	36	8	0.0307	0.0542	855
Hemiculterella wui	24	5	0.0017	0.0643	581
Ischikauia steenackeri	1	1	NA	0.0577	NA
Macrochirichthys macrochirus	1	1	NA	0.1643	NA
Megalobrama amblycephala	28	8	0.0118	0.0057	1076
Megalobrama pellegrini	11	2	0.0000	0	156
Megalobrama skolkovii	8	4	0.0000	0	301
Megalobrama terminalis	48	10	0.0033	0.0299	836
Metzia formosae	2	2	0.0688	0.0368	1624
Metzia lineata	11	5	0.0067	0.0242	477
Metzia longinasus	1	1	NA	0.0367	NA
Metzia mesembrinum	1	1	NA	0.0242	NA
Parabramis pekinensis	47	11	0.0050	0.0471	1656
Pseudohemiculter dispar	57	14	0.1078	0.0913	836
Pseudohemiculter hainanensis	8	3	0.0135	0.0542	688
Pseudolaubuca engraulis	25	11	0.010	0.0892	1259
Pseudolaubuca sinensis	103	15	0.0396	0.0984	1220
Sinibrama macrops	21	6	0.0017	0.0381	962
Sinibrama melrosei	35	11	0.0767	0.0652	748
Sinibrama taeniatus	10	1	0.0151	0.0346	0
Sinibrama wui	1	1	NA	0.0381	NA
Toxabramis houdemeri	74	15	0.0084	0.0646	771
Toxabramis swinhonis	3	1	0.0000	0.0266	0

Mean genetic divergence was 0.93% (0%–6.88%; standard error [SE] = 0.0009) within species and 3.51% (0.67%–5.82%; SE = 0.0034) between species within genera. Maximum intraspecific distances ranged between 0% and 10.78%, and nearest neighbour distances ranged from 0% to 16.43% (Table 1). In addition, eight species displayed maximum intraspecific distances above 2%, including Hemiculter leucisculus (5.56%), Hemiculterella sauvagei (3.07%), Sinibrama melrosei (7.67%), Chanodichthys erythropterus (7.07%), Culter mongolicus (3.96%), Pseudohemiculter dispar (10.78%), Pseudolaubuca sinensis (3.96%) and Metzia formosae (6.88%) (Figure S1; Table 1). Plotting maximum intraspecific and nearest neighbour K2P genetic distances revealed a lack of barcoding gaps with maximum intraspecific genetic distances exceeding minimum distances to the nearest neighbour in multiple cases (Figure 3a). Linear regressions of maximum intraspecific genetic distances as a function of the number of sampled individuals were significant (R = 0.565, p = .001), but had low explanatory power (adjusted R 2 = .058). Maximum intraspecific genetic distances also depended on geographical scale of sampling (adjusted R 2 = 0.143, p = .019) (Figure 3c). Genetic distances to the nearest neighbour declined from an average of 5.10% to 2.81% as the geographical range of sampling was increased from <1 to <5000 km (Figure 3d), with a significant effect (ANOVA; df = 3.934, p = .005).

FIGURE 3

Distribution of genetic distances within and between species and MOTUs, and relationships to geographical scales. (a) Relationship between maximum intraspecific and nearest‐neighbour genetic distance among species. (b) Relationship between maximum intraspecific and nearest‐neighbour genetic distance among MOTUs. Lines indicate the 1:1 ratio. (c) Maximum intraspecific variation against maximum geographical extent (km) of sampled individuals (adjusted R 2 = 0.143, p = .019). (d) the effect of geographical scale of sampling on the closest interspecific divergence. Minimum interspecific divergence across species in five distance categories. In each category, all interspecific distances between individuals with a pairwise geographical distance of less than the category value were calculated. Genetic distance was significantly smaller in the 5000‐km category compared with 1‐ and 50‐km categories (one‐way ANOVA, Tukey's HSD, p < .05)

A visual inspection of the NJ tree based on K2P genetic distances indicated that 18 of 42 species analysed displayed concordant clusters of DNA barcodes (Figure S1). The 24 remaining species encompass multiple lineages, either lumped or mixed due to haplotype sharing. MOTU delimitation analyses yielded varying numbers of MOTUs according to methods with 41, 45, 36 and 53 MOTUs delimited by RESL, ABGD, PTP and GMYC, respectively (Figure 4). MOTU number was underestimated by RESL as BINs were not available for two sequences corresponding to two distinct species in BOLD (Table S1). The final consensus, established by a majority‐rule consensus, consisted of 44 MOTUs. A total of 17 species were unambiguously delimited by the delimitation analyses, eight species displayed multiple MOTUs (Table 2) and 18 species displayed mixed genealogies with nine MOTUs shared by more than one species (Figure 4; Table 3). Distributions of both maximum intraspecific distances and distances to the nearest neighbour separated all observed MOTUs (Figure 3b).

FIGURE 4

DNA‐based species delimitation within the 42 Culterinae species in this study. MOTU delimitation schemes collected from the RESL, ABGD, PTP and GMYC algorithms (yellow), including the consensus delimitation scheme (blue). A Bayesian tree produced in beast version 1.8.2 is used as a base to summarize delimitation results. Scale bar corresponds to Million years

TABLE 2

List of species with multiple MOTUs including their barcode index number in BOLD, maximum genetic distances within an MOTU (D M), and distance to their nearest MOTU (D N). The nearest‐neighbour distance for each MOTU corresponds to the distance to the closest MOTU in this study

	D M	D N
Chanodichthys erythropterus
MOUT4 (BIN:ACU0863)	0.0067	0.0448
MOUT5 (BIN:ADW3458)	0.0000	0.0530
Hemiculter leucisculus
MOUT36 (BIN:ADM2843)	0.0152	0.0390
MOUT39 (BIN:ACB5189)	0.0236	0.0305
Hemiculterella sauvagei
MOUT30 (BIN:ADG7506)	0.0067	0.0242
MOUT31 (BIN:ADC8405)	0.0067	0.0242
Metzia formosae
MOUT42 (BIN:ADC6157)	NA	0.0678
MOUT43 (BIN:ADC6328)	0.0069	0.0678
Pseudohemiculter dispar
MOUT 20 (BIN:AEE2041)	0.0084	0.0897
MOUT 21 (BIN:AEH4884)	0.0084	0.0521
MOUT 22 (BIN:AEH7332)	0.0084	0.0434
MOUT 23 (BIN:AEH7333)	NA	0.0434
Pseudohemiculter hainanensis
MOUT33 (BIN:AEH4806)	0.0000	0.0121
MOUT34 (BIN:AEH4807)	0.0017	0.0121
Pseudolaubuca sinensis
MOUT25 (BIN:ADC6253)	0.0203	0.0343
MOUT26 (BIN:ACU0850)	0.0033	0.0343
Sinibrama melrosei
MOUT15(BIN:AEH7279)	0.0050	0.0571
MOUT16 (BIN:AEH3997)	NA	0.0567
MOUT17 (BIN:AEH3998)	0.0101	0.0335
MOUT18 (BIN:AEH8015)	0.0050	0.0335

TABLE 3

List of MOTUs shared by multiple species, including their maximum genetic distances within an MOTU (D M), distance to their nearest MOTU (D N), and a list of species detected for each MOTU. The nearest‐neighbour distance for each MOTU corresponds to the distance to the closest MOTU without considering species boundaries

	D M	D N	Species
MOTU1	0.0203	0.0315	Culter alburnus, Culter mongolicus, Cuter recurviceps, Chanodichthys erythropterus
MOTU3	0.0033	0.0328	Ancherythroculter nigrocauda, Culter dabryi,
MOTU4	0.0067	0.0448	Chanodichthys compressocorpus, Chanodichthys erythropterus
MOTU7	0.0017	0.0233	Ancherythroculter kurematsu, Ancherythroculter wangi
MOTU10	0.0118	0.0314	Megalobrama Amblycephala, Megalobrama pellegrini, Megalobrama skolkovii
MOTU29	0.0101	0.0480	Anabarilius liui chenghaiensis, Anabarilius liui yalongensi
MOTU35	0.0117	0.0602	Hemiculter bleekeri, Hemiculter lucidus
MOTU39	0.0236	0.0305	Hemiculter leucisculus, Toxabramis swinhonis
MOTU43	0.0069	0.0678	Metzia formosae, Metzia longinasus

Haplotype networks were reconstructed for 18 species displaying shallow genetic divergence and/or haplotype sharing (Figure 5). Five cases of haplotype sharing were observed between Culter alburnus and Cu. mongolicus (Figure 5a), Ancherythroculter nigrocauda and Culter dabryi (Figure 5b), Chanodichthys erythropterus and Ch. compressocorpus (Figure 5d), Hemiculter bleekeri and H. lucidus (Figure 5e), and Megalobrama skolkovii and M. pellegrini (Figure 5g). For each of these five cases, haplotype sharing involves ancestral haplotypes, as revealed by their central position in haplotype network reconstructions. Scattered haplotypes across haplotype networks were observed for Ch. erythropterus (Figure 5a), with haplotypes originating from Hainan Island being more closely related to haplotypes of Culter recurviceps than conspecific haplotypes sampled in the drainages of China's mainland.

FIGURE 5

Median‐joining networks of DNA barcode haplotypes for the species groups involved in mixed genealogies and/or haplotype sharing. (a) Three Culter species and Chanodichthys erythropterus. (b) Ancherythroculter nigrocauda and Culter dabryi. (c) Ancherythroculter kurematsui and A. wangi. (d) Ch. Erythropterus and Ch. Compressocorpus. (e) Hemiculter bleekeri and Hemuculter lucidus. (f) Hemiculter leucisculus and Toxabramis swinhonis. (g) Three Megalobrama species. (h) Metzia formosae and M. longinasus. Empty circles and numbers indicate unsampled haplotypes and mutation steps, respectively

NJ trees built for species displaying multiple, and deeply diverging, MOTUs revealed a diversity of patterns. Two lineages were detected in Hemiculter leucisculus, Ch. erythropterus, Pseudolaubuca sinensis, Hemiculturella sauvagei, Pseudohemiculter hainanensis and Metzia formosae and four lineages were observed in Pseudohemiculter dispar and Sinibrama melrosei (Figure 6). In most cases of two MOTUs, a south–north differentiation was observed, as exemplified by Ch. erythropterus, Hemiculter leucisculus, Hemiculturella sauvagei and Pseudolaubuca sinensis, with one lineage distributed in the southern drainages including Hainan Island, Pearl River, and some coastal rivers in southern China, and another lineage largely occurring in central and northern drainages such as the Yangtze River and/or its northern drainages (Figure 6). However, the co‐occurrence of haplotypes from distinct MOTUs is observed in some tributaries of the Yangtze and the Pearl River (Figure 6a–d). Alternatively, allopatric distributions of conspecific MOTUs involved different patterns with Vietnam vs. Taiwan distribution of MOTUs for M. formosae, continental vs. Hainan Island MOTUs for Pseudohemiculter hainanensis, and erratic distribution in the case of Pseudohemiculter dispar and S. melrosei (Figure 6). Divergence time estimates suggest most MOTU divergence events originated during the Pleistocene, and several exceptions are detected within Pseudohemiculter dispar, S. melrosei and M. formosae with MOTU divergence dated back to the late Pliocene (Figure S2).

FIGURE 6

Phylogeographical patterns of eight species with multiple lineages and/or MOTUs. (a) Chanodichthys erythropterus. (b) Hemiculter leucisculus. (c) Hemiculterella sauvagei. (d) Pseudolaubuca sinensis. (e) Metzia formosae. (f) Pseudohemiculter hainanensis. (g) Pseudohemiculter dispar. (h) Sinibrama melrosei. Trees at the top left of each map are neighbour‐joining trees of the corresponding species K2P model. Different colour circles represent different lineages and/or MOTUs. Scale bars correspond to K2P genetic distances

DISCUSSION

The present study provides the first comprehensive DNA barcode reference library of the subfamily Culterinae, including 42 species (~66% of known Culterinae species) and covering most major East Asian watersheds (Luo & Chen, 1998). Notably, all Anabarilius species and several monotypic genera such as Pogobrama and Hainania are missing. These species have very restricted distribution ranges and little information is available regarding their abundance and ecology (Jiang et al., 2016; Luo & Chen, 1998). In addition, DNA sequences are scarce for these species in international repositories such as GenBank, EMBL and BOLD. By contrast, most widespread species were aptly covered, with multiple individual sequences originating from different East Asian drainages, which allows assessment of the impact of geographical scale on the performance of DNA barcoding in this subfamily.

Shallow divergence, haplotype sharing and taxonomic confusion

DNA‐based species delimitation analyses converged with morphological delimitation in only 17 species, indicating COI sequences fail at capturing species boundaries in multiple Culterinae. This relatively low success rate (~40%) was unexpected, considering the large number of empirical studies which yielded better matches between morphological and DNA‐based delimitations (>90%) than here (April et al., 2011; Dahruddin et al., 2017; Delrieu‐Trottin et al., 2019; Hou et al., 2018; Hubert et al., 2008; Keskin & Atar, 2013; Lara et al., 2010; Pereira et al., 2013; Ribeiro et al., 2012; Shen et al., 2016; Sholihah et al., 2020; Sonet et al., 2019; Ward et al., 2005). The lower success rate observed here is probably the result of the broad spatial scale involved and high coverage of closely related species in poorly known taxa. This is exemplified by the lower mean interspecific distance (3.51%) than observed in previous studies focused on fish families or subfamilies within narrow geographical areas (Adeoba et al., 2018; Alcantara & Yambot, 2016; Hubert et al., 2010; Jaafar et al., 2012; Tzeng & Chiu, 2012; Zhang et al., 2017). This trend is further supported by the statistically significant effect of geographical distance on both maximum intraspecific genetic distance (positively correlated) and genetic distance to the nearest neighbour (negatively correlated) (Figure 3). This result suggests the presence of a substantial proportion of recently diverged species among Culterinae species (Bellafronte et al., 2013), a trend further supported by the discovery of 18 cases of mixed interspecific genealogies (Culter spp., Chanodichthys spp., Hemiculter spp., Megalobrama spp., Metzia spp. and Toxabramis swinhonis). However, mismatches between morphological and DNA‐based species delimitations may arise from introgressive hybridization and incomplete lineage sorting, but taxonomic confusion and erroneous identifications cannot be discarded (Dahruddin et al., 2017; Durand et al., 2017; Hubert & Hanner, 2015; Sholihah et al., 2020; Sonet et al., 2019).

In the genus Megalobrama, three species (i.e., M. amblycephala, M. skolkovii, and M. pellegrini) are molecularly indistinguishable, as they share the same mitochondrial lineage (MOTU10, Figure 3), and haplotype sharing was uncovered between M. skolkovii and M. pellegrini (Figure 5g). These close genetic relationships were previously mentioned (Bai et al., 2015; Xie et al., 2012), but no haplotype sharing was observed at the mitochondrial ND2 gene (Bai et al., 2015), suggesting this lack of resolution results from the lower substitution rate at the COI than ND2 gene. Interspecific hybridization experiments demonstrated that Megalobrama species can hybridize (Zhang et al., 2014), but M. skolkovii and M. pellegrini have allopatric distributions. M. skolkovii occurs in the middle and lower Yangtze River, and the Pearl River, and M. pellegrini occurs in the upper Yangtze River (Luo & Chen, 1998). A similar pattern is observed within the genera Culter, Hemiculter, Ancherythroculter and Chanodichthys with haplotype sharing occurring between species pairs (Cu. mongolicus and Cu. recurviceps, Figure 5a; Ancherythroculter nigricauda and Culter dabryi, Figure 5b; Chanodichthys erythropterus and Ch. compressocorpus, Figure 5d; Hemiculter bleekeri and Hemuculter lucidus, Figure 5e). In these particular cases, shared haplotypes have a central position in reconstructed haplotype networks, suggesting these are ancestral haplotypes. In addition, no haplotype sharing is detected among recently diverged haplotypes, suggesting that recent divergence of mitochondrial lineages is more probably responsible for this shortcoming than introgressive hybridization (Hubert & Hanner, 2015). Along the same line, taxonomic gaps or erroneous identifications are unlikely because of the presence of diagnostic morphological characters such as mouth shape, ventral ridge, chamber number and body type for most of the species pairs involved in haplotype sharing (Luo & Chen, 1998). This is particularly evident in Ancherythroculter nigrocauda and Cu. dabryi, which have been assigned to distinct genera due to different numbers of chamber in the swim bladder, Ancherythroculter spp. and Culter spp. having two and three chambers, respectively (Ding, 1994; Luo & Chen, 1998; Zheng, 1989). By contrast, Hemiculter bleekeri and Hemuculter lucidus are distinguished by quantitative characters including head–body length ratio (Luo & Chen, 1998). Moreover, Hemuculter lucidus was previously considered as a subspecies of Hemiculter bleekeri due to their morphological similarity (Luo & Chen, 1998; Wu & Yi, 1959). Considering that there is no overlap in the distribution of the two species (Luo & Chen, 1998), taxonomic confusion is more likely in this case and the two species should probably be regarded as a single species.

Cryptic diversity and phylogeographical patterns

The mean level of intraspecific divergence observed here (0.93%) is higher than most of those previously reported in DNA barcoding studies focusing on fish families or subfamilies (0.75%) (Adeoba et al., 2018; Hubert et al., 2010; Jaafar et al., 2012; Laskar et al., 2013; Shen et al., 2016; Tzeng & Chiu, 2012; Zhang et al., 2017). The high intraspecific divergence observed in Culterinae is probably explained by the large geographical scale of the present study, widespread species being mostly spatially structured (Avise, 2000; Bergsten et al., 2012; Nekola & White, 1999; Wright, 1943), a trend previously observed in other Cypriniformes subfamilies such as Rasborinae in Southeast Asia (Sholihah et al., 2020). Several cases of multiple, highly divergent intraspecific MOTUs were detected in eight species, including Chanodichthys erythropterus, Hemiculter leucisculus, Hemiculterella sauvagei, Metzia formosae, Pseudolaubuca sinensis, Pseudohemiculter spp. and Sinibrama melrosei (Figure 4). Most intraspecific lineage divergences are dated to between the Pliocene and Pleistocene (Figure S2), suggesting the influence of historical geological events.

In Ch. erythropterus, Hemiculter leucisculus, Hemiculturella sauvagei and Pseudolaubuca sinensis, high divergence (>3.2% maximum intraspecific genetic distance) between MOTUs is observed. For Ch. erythropterus and Hemiculter leucisculus, MOTUs display alternative geographical distributions in the southern coastal watersheds, or the Pearl River and the Yangtze River. Similar distribution patterns were previously reported in phylogeographical studies of Hemiculter leucisculus and Opsariichthys bidens (Chen et al., 2017; Perdices et al., 2005). For Hemiculturella sauvagei and Pseudolaubuca sinensis, two MOTUs were identified from the Pearl River and Yangtze River sampling sites, with one MOTU in the Yangtze and the other in both the Pearl River and southern tributaries of the Yangtze River (Figure 6). These distribution patterns are consistent with the influence of the Nanling Mountains, elevation of which created a biogeographical barrier separating the Yangtze River and the Pearl River, a scenario previously proposed by phylogeographical studies of multiple fish taxa, such as Hemiculter leucisculus (Chen et al., 2017), O. bidens (Li et al., 2009; Perdices et al., 2005), Glyptothorax fokiensis (Chen et al., 2007) and Siniperca scherzeri (Wang et al., 2006). Elevation of the Nanling Mountains occurred during uplift of the Qinghai–Tibetan plateau between 3.6 and 0.15 million years ago (Ma), a time frame matching our divergence time estimates among these MOTUs, dated at 1.36–2.32 Ma (Figure S2). Interactions between geology and palaeoclimatic fluctuations are also likely, as the Asian monsoon systems became established at 3.6–2.6 Ma. These climatic and geological events probably drove landscape features and drainage systems of eastern Asia, promoting geographical isolation and subsequent divergence in many fish species (Chen et al., 2007; Chen et al., 2017; Guo et al., 2006; He et al., 2001; Perdices et al., 2005; Wang, 2005; Yang et al., 2009; Yu et al., 2014). In addition, several sequences of Hemiculturella sauvagei and Pseudolaubuca sinensis in some southern tributaries of the Yangtze River were grouped with the Pearl River lineage, and many sequences of Hemiculter leucisculus and Ch. erythropterus from the northern tributaries of the Pearl River were clustered with sequences of the Yangtze River and its northern tributaries. This is a common pattern observed in many fish taxa distributed in neighbouring tributaries of these two river basins (Chen et al., 2007; Chen et al., 2017; Perdices et al., 2005; Wang, 2005; Yang et al., 2009) due to intricate landscape changes involving watershed reconfiguration through headwater capture events, which may have thereby promoted fish dispersal by temporarily connecting upper branches of the Pearl and Yangtze Rivers (Ren et al., 1959).

With regard to Pseudohemiculter dispar and Sinibrama melrosei, four observed lineages were detected and did not show marked geographical distribution patterns (Figure 6g,h), suggesting these two species may have initially undergone allopatric divergence, followed by secondary contacts due to drainage rearrangements occurring in southern China. The above‐mentioned climatic and geological events during the Pliocene and Pleistocene may have contributed to generate multiple MOTUs. In addition, previous studies demonstrated that geomorphological changes of drainage networks in response to the late Pleistocene glaciations had occurred in southern China (Liu et al., 2003; Wang et al., 1999), providing many opportunities for fish dispersal.

A repeated pattern of alternative distribution of MOTUs or haplotypes was also observed between Hainan Island and the mainland in Pseudohemiculter hainanensis (Figure 6f). This pattern suggests the formation of the Qiongzhou Strait triggered the divergence of populations within Pseudohemiculter hainanensis during the mid‐Pleistocene (~0.53 Ma). The influence of this strait on population divergence has been previously suggested for several freshwater fish species such as Garra orientalis (Yang et al., 2016), Micronemacheilus pulcher (Qiu et al., 2008), Megalobrama terminalis (Chen et al., 2020) and Glyptothorax hainanensis (Chen et al., 2007).

Implications for Culterinae taxonomy and future prospects

Our study provides the most comprehensive geographical, taxonomic and molecular sampling of this subfamily to date, with several findings contradicting current systematic and taxonomic knowledge. The most challenging cases were those involving shallow interspecific divergence between species, a trend shared by 18 species. Given that most of these cases may arise from recent divergence during the Pleistocene, additional molecular markers with higher substitution rates, such as the mitochondrial ND2 gene and control region, would certainly help in better distinguishing cases of recent divergence from actual haplotype sharing. Nevertheless, haplotype sharing was also observed between Ancherythroculter nigrocauda and Culter dabryi, for instance, calling for a broader assessment of the molecular species boundaries using biparentally inherited nuclear loci, and resulting in the revision of the generic status of the two species. With regard to Hemiculter bleekeri and Hemuculter lucidus, the relevance of the head–body length ratio in distinguishing each species should be reconsidered. However, we cannot rule out the occurrence of introgressive hybridization between these closely related species. Introgression is known to distort the congruence between morphological and DNA‐based species delimitations (Petit & Excoffier, 2009), and constitutes a challenge for species identification based on DNA barcodes. Using additional nuclear loci in combination with COI in a spatially explicit context will probably shed light on species identity and ongoing evolutionary dynamics driving their genetic divergence (Oliveira et al., 2010; Raupach et al., 2010).

At the other end of the spectrum, the detection of multiple, and highly divergent, MOTUs in eight species suggests Culterinae diversity is currently underestimated and potentially new Culterinae species are awaiting a formal description. This pattern can be explained by overlooked, subtle morphological differences, but cryptic diversity and unrecognized speciation events might also be responsible for this result (Hebert et al., 2004; Hubert et al., 2012). The taxonomic status of these cryptic lineages delimited by mitochondrial sequences is unstable (Satler et al., 2013). These putative species constitute primary species hypotheses awaiting further assessments of their morphological divergence and phylogenetic placement using additional nuclear loci to produce secondary species hypotheses and robust species descriptions.

CONCLUSION

The present large‐scale genetic study in East Asia demonstrates the limit of DNA barcoding for automated identification of Culterinae species, as DNA barcodes and morphological characters are congruent for only 17 species. The 18 cases belonging to a complex of closely related species currently challenge the implementation of routine molecular species identification, but also question the evolutionary origin of a substantial number of Culterinae species. Moreover, newly detected mitochondrial lineages in nine species suggest Culterinae diversity might be underestimated, which calls for a substantial effort to better characterize the diversity of this particular group and has important implications in terms of conservation and management. This result further indicates East Asia provided important evolutionary reservoirs during the diversification of Culterinae. As such, this study warrants further research developments and provides guidelines for future taxonomic and genomic studies.

AUTHORS CONTRIBUTIONS

W.C., N.H., X.L. and J.L. designed the study; W.C., D.X. and C.Z. performed molecular experiments; W.C., Y.L., X.C., S.Z., J.Y., X.W. and C.Z. collected specimens and coordinated morphological identifications. W.C. and N.H. performed data analyses and wrote an initial draft of the manuscript. All authors contributed to the writing of the final draft.

CONFLICT OF INTEREST

The authors declare no competing interests.

OPEN RESEARCH BADGES

This article has obtained an Open Data Badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5883/DS‐CUBAR.

Figure S1

Click here for additional data file.

Figure S2

Click here for additional data file.

Table S1

Click here for additional data file.