Phylogeographic analysis revealed allopatric distribution pattern and biogeographic processes of the widespread pale chub Opsariichthys acutipinnis-evolans complex (Teleostei: Cyprinidae) in southeastern China

1. Introduction

Understanding the spatial arrangements of genetic variations within populations or among closely related species is a major objective of phylogeography (e.g., Avise et al., 1987; Avise, 2009). Interpretations of phylogeographic patterns and genetic structures provide insights into speciation processes and guidance for conservation and management decisions (Avise, 2000). For example, population genetic and phylogeographic studies revealed that the widespread but endangered Chinese giant salamander (Andrias davidianus) actually harbors seven distinct genetic lineages that might represent at least five different species (Murphy et al., 2000; Yan F. et al., 2018; Liang et al., 2019; Chai et al., 2022), and the lack of awareness of the existence of these cryptic lineages misguided conservation efforts in the past decades. In contrast, genetic assessment of the economically important Lake Tanganyika sprat Stolothrissa tanganicae indicates this species lacks population structure and suggests an integrated stock management strategy (De Keyzer et al., 2019). Understanding the genetic variations and structures of many other species also enhanced our knowledge of the divergence patterns and offered management suggestions (e.g., Van Rossum et al., 2018; Guo et al., 2019; Rincon-Sandoval et al., 2019; Tumendemberel et al., 2019). These findings further emphasize the necessity of genetic evaluations to determine the taxonomic status and identify evolutionary significant units (ESUs) in conservation initiatives, especially for seemingly widespread species (e.g., Taylor et al., 2011; Larson et al., 2014; Yan F. et al., 2018; Li et al., 2020; Li C. et al., 2023).

Opsariichthys acutipinnis and O. evolans are two small riverine fishes widely distributed throughout southeastern China, including the Yangtze River basin, Pearl River basin and independent coastal river basins like the Qiantangjiang River, Minjiang River, and Hanjiang River. These fishes prefer lotic environments with gravel and pebble substrates, and are most common in mountain streams and small tributaries (Chen and Chu, 1998). Opsariichthys acutipinnis was first described by Bleeker in 1871 as Barilius acutipinnis based on a specimen from the Yangtze River (Bleeker, 1871), and then being transferred to the genus Opsariichthys without any explanation (Günther, 1873). Later in 1902, Jordan and Evermann established the genus Zacco to include fishes originally assigned to the genus Opsariichthys, but without conspicuous notched jaws (Jordan and Evermann, 1902). They also described Zacco evolans from the Taiwan Island. Later studies generally treated Opsariichthys acutipinnis and Opsariichthys evolans as synonyms of Opsariichthys platypus (Rendahl, 1928) or Zacco platypus (Oshima, 1919; Nichols, 1943; Yang and Huang, 1964), but did not reach any consensus. In 1982, Chen (1982) revised the taxonomy of opsariichthyine fishes and agreed that the notched jaws could distinguish the genus Opsariichthys and Zacco, thus he considered both O. acutipinnis and O. evolans should be member of the genus Zacco and be synonymized with Zacco platypus. However, some later studies based on osteological or molecular evidence suggest that the notched jaws might not be appropriate to serve as diagnostic feature (Howes, 1980; Ashiwa and Hosoya, 1998; Chen and Chang, 2005; Ma et al., 2006). Chen et al. (2009) suggested that the genus Zacco can be differed from Opsariichthys by vertical bars on flanks fused into patches (versus distinct vertical bars not fused into patches in Opsariichthys) and nuptial tubercles on cheeks in male fused into a plate basally (versus tubercles small and not fused into a plate in Opsariichthys). Thus the species previously identified as Zacco platypus from southern China with such features should be regarded as members of Opsariichthys acutipinnis-O. evolans complex. Later studies generally treated O. acutipinnis and O. evolans as two closely related but distinct species based on deep divergence in mitochondrial gene (s), allopatric distribution, and morphological differences (Huynh and Chen, 2013; Wang et al., 2019). A recent work conducted by Wang et al. (2019) revealed a deeply divergent mitochondrial lineage that related to O. acutipinnis and O. evolans (“O. sp. A”), indicating the existence of previously ignored cryptic diversity. In the past decade our extensive sampling work throughout southeastern China has yielded plenty of morphologically similar specimens that could be identified as either O. acutipinnis and O. evolans. However, we found out that morphological identification based on diagnostic characteristics proposed by Huynh and Chen (2013) and Wang et al. (2019) would be problematic, e.g., specimens with features of O. acutipinnis and O. evolans co-existed in most populations. As a consequence, we tentatively treated these fishes as the Opsariichthys acutipinnis-evolans complex in the current study.

Plenty of studies have suggested that primary freshwater fishes with resident lifestyles always exhibit highly structured genetic patterns between different drainage basins or even different stretches of a river because of their typically low mobility and preference for certain habitats (e.g., Chen et al., 2007; Lin et al., 2008; Huang et al., 2019; Won et al., 2020). For example, populations of the cold-water fish Rhynchocypris oxycephalus in China exhibited deep genetic structure, aside from the geological factors, the low dispersal ability and preference for headwater streams are important contributors to the current divergence patterns (Yu et al., 2014). Prominent genetic differentiation was detected from the stream goby populations within fine-scale regions, and that both paleodrainage alternations and the sedentary lifestyle of this species are major factors in generating the strong genetic structure (Wu et al., 2016). Cases like these reinforce the impression that primary freshwater fishes are ideal candidates to investigate the biological consequences of complex geological changes. The O. acutipinnis-evolans complex is a group of small riverine fishes widely distributed throughout southeastern China, therefore they can be great models to explore how geological and climatic changes affect their genetic makeups and distribution patterns. However, few studies (Perdices et al., 2004; Perdices and Coelho, 2006; Lin, 2008; Yin, 2015) have addressed the genetic differentiation and phylogeographic structure of this species complex due to long term neglect in taxonomy, which hindered us from understanding the diversity and evolutionary history of these small stream fishes, as well as making proper management decisions. In the present study, we used both mitochondrial Cyt b sequences and 25 morphological traits to assess the genetic and phenotypic differentiation among populations of O. acutipinnis-evolans complex in southeastern China. Then we aim to investigate the phylogeographic structure and explore the potential factors that have shaped the patterns of genetic variability. Finally, we wish our phylogeographic evaluation could offer some implications for the conservation and management of this species.

2. Materials and methods

2.1. Sampling, DNA extraction, PCR amplification, and sequencing

A total of 821 specimens were collected from 36 localities in southeastern China, covering the Yangtze River basin, Pearl River basin, and several small coastal river basins such as the Qiantangjiang River, the Minjiang River and the Hanjiang River (Supplementary Table S1). Fishes were captured using dip net, gill net, and electrofishing between 2009 and 2020. Due to the lack of a generally acceptable classification system, we treated all our Opsariichthys specimens with the following characteristics as Opsariichthys acutipinnis-evolans complex: (1) mouth without conspicuous notch; (2) lateral line scale count lower than 45; (3) tubercles round and small, usually arranged in one row on the lower jaw. Voucher specimens were preserved in 95% ethanol at 4°C. All voucher specimens were deposited at the Institute of Hydrobiology, Chinese Academy of Sciences.

Total genomic DNA was extracted either from a small piece of flesh or pelvic fin clip taken from the right side of each specimen, following the salt-extraction methods (Aljanabi, 1997). The complete mitochondrial cytochrome b (Cyt b) gene was amplified for all the samples with the universal primers L14724 and H15915 (Xiao et al., 2001). Amplification was conducted in a 30-μL reaction volume, containing 11 μl of dd H₂O, 15 μl of 2 × MIX, 1 μl of each primer, and 2 μl of template DNA. The thermocycling program for the polymerase chain reaction (PCR) consisted of an initial denaturation at 94°C for 5 min, followed by 35 cycles of denaturation at 94°C for 45 s, annealing at 54°C for 45 s, extension at 72°C for 1 min, and a final extension at 72°C for 10 min. PCR products were sent to commercial sequencing company (Sangon Biotech, Shanghai) for purifying and sequencing with the same PCR primers.

We sequenced 622 samples and retrieved 199 sequences from our previous study (Wang et al., 2019). We also included another 193 sequences from GenBank submitted by previous authors (Perdices et al., 2004; Perdices and Coelho, 2006; Lin, 2008; Zheng et al., 2015; Yin et al., 2016; Wang et al., 2019; Supplementary Table S1). Thus a total of 1014 Cyt b sequences from 62 localities were used in the present study, covering most of its distribution in southeastern China, including the Yangtze River basin, Pearl River basin, and several small coastal river basins (Figure 1A, Table 1, Supplementary Table S1). As the information on the taxonomic status of some sequences from previous studies have not been updated, those downloaded from GenBank were verified using sequences from Wang et al. (2019) as reference to confirm their identity prior to the analysis. Opsariichthys macrolepis (OQ428185, OQ428186), O. hainanensis (OQ428183, OQ428184), and O. chengtui (KT725244)were used as outgroup taxa. All the newly determined sequences were deposited in GenBank under the accession numbers OQ214215 to OQ214836.

Figure 1. (A) Map of sampling localities of the O. acutipinnis-evolans complex. Codes for sampling localities are given in Table 1 and Supplementary Table S1. Populations are present as pie charts with slice size proportional to the frequency of the six major lineages identified by the Cyt b gene tree as in (B). (B) Cyt b gene tree from Bayesian inference. Numbers in bold and italics at nodes indicate Bayesian posterior probabilities (PP) and bootstrap values (BS) from maximum likelihood analysis, respectively. Only well-supported numbers (PP > 0.95, BS > 70) are shown. Numbers in regular at nodes indicate divergence time estimate of major lineages. Colored branches denote different lineages as in (A). (C–F) Median-joining (MJ) networks of the O. acutipinnis-evolans complex concerning lineage (B, D-F). The mutation steps are shown as hatch marks on the branches. Missing (or unsampled) haplotypes are represented by black dots. The area of each circle is proportional to the number of individuals for each haplotypes. Different colors denote different populations.

Table 1. Sampling information for the 62 populations of O. acutipinnis-evolans complex used in this study.

2.2. Sequence analyses

Multiple sequence alignment was performed using MAFFT (Katoh and Standley, 2013) implemented in PhyloSuite (Zhang et al., 2020) under default parameters and subsequently checked by eye in SeaView (Gouy et al., 2010). Characteristics and variations of sequences were analyzed in MEGA 7.0 (Kumar et al., 2016). Unique haplotypes were collapsed using ALTER (Glez-Peña et al., 2010). Sequence diversity indices were estimated for both pooled and individual populations using Arlequin 3.5 (Excoffier and Lischer, 2010), including the number of haplotypes (n), haplotype diversity (h) and nucleotide diversity (π).

2.3. Phylogenetic analyses

Phylogenetic relationships among mitochondrial Cyt b haplotypes were reconstructed using Bayesian inference (BI) and maximum likelihood (ML) approaches. TN + F + G4 was selected as the best-fit nucleotide substitution model based on the corrected Akaike information criterion (AICc) in Modelfinder (Kalyaanamoorthy et al., 2017).

BI was conducted in MrBayes 3.2.7 (Ronquist et al., 2012) with two independent runs of Metropolis-coupled Markov Chain Monte Carlo analyses, each comprising one cold chain and three heated chains at a default temperature of 0.1. The chains were run for 10 million generations and sampled every 1,000 generations. The average standard deviation of split frequencies were required to drop below 0.01, and the convergence diagnostic for branch length posterior probabilities (potential scale reduction factor) to approach 1. The first 25% of the trees were discarded as burn-in and the remaining trees were used to generate a consensus tree. ML analysis was performed in PhyML 3.0 (Guindon et al., 2010). The robustness of each branch was estimated with 1000 non-parametric bootstrap replicates. To visualize the relationships among haplotypes in each lineages, the median-joining (MJ) networks were generated using PopArt (Leigh and Bryant, 2015) on lineages that contained multiple populations.

2.4. Estimation of divergence time

A dated gene tree using the Cyt b haplotypes was reconstructed using BEAST 2.6.0 (Bouckaert et al., 2019) to estimate divergence time among major lineages. The hypothesis of a strict molecular clock was first tested using the sequence-based likelihood ratio test in DAMBE 7(Xia, 2018). The hypothesis of a strict molecular clock cannot be rejected at the significance level of 0.05 (p = 0.8798). We calibrated the Cyt b gene tree with 1.05% per site per million year under the strict clock model to obtain the absolute divergence time estimates. This substitution rate was previously inferred for several opsariichthyine fishes from Lake Biwa based on geological event calibration (Tabata et al., 2016). Analyses were performed for 10 million generations while sampling every 1000 generations under the Yule speciation tree prior. The BI tree was used as the starting tree. Three independent runs were carried out. Tracer 1.7 (Rambaut et al., 2018) was used to check for adequate effective sample sizes (ESS > 200) of the posterior distribution of each parameter and to access convergence among independent runs. The runs were then combined using LogCombiner (Bouckaert et al., 2019) after removing the first 25% samples as burn-in. A maximum clade credibility tree (MCC tree) was generated using TreeAnnotator in the BEAST package.

2.5. Morphological examination

A total of 273 specimens representing each lineages were randomly selected for morphological examination. Twenty morphometric measurements and three meristic characters were used in this study (Supplementary Table S2). Besides, the position of the posterior end of pectoral fin and maxilla were check visually. Measurements and counts were made on the left side of each individual following Armbruster (2012). Measurements were taken point to point with a digital caliper (precision of 0.01mm). The pectoral fin was removed from the selected specimen, cleared and stained following Taylor and van Dyke (1985), and then inspected under a stereoscope.

Morphometric measurements were subject to principal component analysis (PCA) to check body shape differentiations among major lineages. Prior to the analysis, the method of Reist (1985) was used to normalize all the measurements except for standard length to eliminate the influence of allometry of body parts. Normalization was conducted in R and PCA was run with SPSS 21 (SPSS, Chicago, IL, USA). The PCA loadings were presented in Supplementary Table S3.

3. Results

3.1. Sequence characteristics

The final alignment of Cyt b consisted of 912 base pairs (bp), of which 174 bp (19.1%) were parsimony informative. The base composition was unequal (A = 24.29%, T = 30.75%, G = 16.44%, C = 28.52%), with strong compositional biases against G existed at the second and especially at the third codon position (just an average of 10.5%). We translated all sequences into amino acids, and no internal stop codon was detected. A total of 263 haplotypes were identified among the 1014 Cyt b sequences (Table 1), of which 84.4% (222) were unique to certain populations. The Cyt b haplotype diversity (h) and nucleotide diversity (π) varied among populations (Table 1). The overall h in all samples was 0.9655 ± 0.0028, and the overall π was 0.157542 ± 0.076141.

3.2. Cyt b gene tree and haplotype networks

The ML and BI trees based on Cyt b haplotypes yielded similar gene tree topologies and congruently revealed six evolutionary lineages clearly correlating with their geographical distribution (Figure 1B). Lineage A located at the basal part of the tree, comprised all samples from population 6–8 (Tt, Xj and DQ) belonging to the Lingjiang River basin on China’s southeastern coast. Lineage B, sister to the remaining four lineages, consisted of samples from population 40–62 belonging to the middle Yangtze River basin and the Pearl River basin. Lineage C consisted of samples from population 39 (JD) of the Qingyijiang River belonging to the lower Yangtze. Lineage D contained samples from northern Taiwan island (population 36–37) and northern populations in China’s southeastern coast (population 2–5, 9–12) except for one population (population 1, AJ), which originated from Tiaoxi River of the lower Yangtze River. Lineage E and F contained samples from middle and southern populations (population 13–35) in the southeastern coastal region, respectively. Population 38 (Dh) from western Taiwan island clustered within lineage F. The monophyly and relationships of lineages was well supported (PP>0.95, BS>75) unless mentioned otherwise. The relationships among haplotypes in lineage B, D, E, and F were clearly shown by the MJ networks (Figures 1C-F). It was clearly depicted that the northern Taiwan populations (population 36, Ss, and population 37, Qc) were derived from Qiantangjiang/Oujiang River poopulations, and the northwestern Taiwan population (population 38, Dh) originated from Mulanxi River population. The inter-lineage genetic distances ranged from 0.0177 (between lineage E and F) to 0.0700 (between lineage A and C) (Table 2).

Table 2. Genetic distance between major lineages of O. acutipinnis-evolans complex based on Cyt b sequences.

3.3. Divergence time estimation

The divergence time estimates of major lineages are provided in Figure 1B. The dating analysis suggested that the most recent common ancestor (TMRCA) of O. acutipinnis-evolans complex was dated back to 3.78 Ma (95% HPD 3.10 ~ 4.44 Ma). TMRCA of lineage A occurred at 3.00 Ma (95% HPD 2.44 ~ 3.56 Ma). TMRCA of lineage B was dated to 2.61 Ma (95% HPD 2.12 ~ 3.11 Ma). TMRCA of lineage C diverged at 2.12 Ma (95% HPD 1.67 ~ 2.59 Ma), and TMRCA of lineage D diverged at 1.33 Ma (95% HPD 1.02–1.66 Ma). The split between the continental and northern Taiwan island populations occurred at 0.84 Ma (95% HPD 0.61 ~ 1.10 Ma). Time since the split of lineage E and lineage F was estimated to be 0.95 Ma (95% HPD 0.70 ~ 1.21 Ma). The northwestern Taiwan population was isolated from the continental populations at around 0.24 Ma (95% HPD 0.13 ~ 0.35 Ma).

3.4. Morphological examination

The scatter plot of the first two principal components indicated no clear separation between the six evolutionary lineages (Figure 2A). The number of scales and branched pectoral fin rays counts showed subtle differences and significant overlap among different lineages (Figures 2B–D; Supplementary Table S2). The most frequently observed predorsal scales ranged from 15 to 16, and circumpeduncular scales ranged from 16 to 18. The clear and stain procedure revealed that the branched pectoral fin rays ranged from 13 to 16, with 14 and 15 being the most frequent counts. Besides, the pectoral fin length varied among individuals, sex and life stages. Longer pectoral fins reaching or extending beyond the origin of the pelvic fin were only observed in sexually active males. The posterior end of pectoral fin never extends to pelvic fin origin in females and juveniles (Figures 2E,F; Supplementary Table S2). The position of the posterior end of maxilla also varied among individuals and did not exhibit specific patterns between lineages (Figure 2G; Supplementary Table S2).

Figure 2. (A) Scatter plot of the first two principal components based on 20 measurements. The meristic counts differences between different lineages in (B) predorsal scales, (C) circumpeduncular scales, and (D) branched pectoral fin rays. The differences in the position of posterior end of pectoral fin against ventral fin origin between different (E) lineages and (F) sexes. (G) The differences in the position of the posterior end of maxilla against vertical of anterior margin of orbit. N, not reaching; R, reaching; EB, extend beyond.

4. Discussion

4.1. Biogeographic processes of Opsariichthys acutipinnis-evolans complex in southeastern China

4.1.1. Geological events drove the early diversification

One of the major goals of phylogeography is to decipher the spatial arrangements of genetic variations within populations or among closely related species (Avise, 2000). Our results clearly revealed that phylogeographic pattern of the O. acutipinnis-evolans complex in southeastern China corresponded to category I defined by Avise (2000), which was allopatric lineages with deep differentiation. This pattern could be explained by the existence of long-term biogeographic barriers. While lineage A, C, D, E, and F originated from the relatively eastern part of their range, including the coastal area of southeastern China, lineage B occupied the relatively western part, covering the middle Yangtze River basin and the Pearl River basin. The most conspicuous geographic break between the eastern and western lineages is the Huangshan Mountains and the Wuyi mountains. These mountains serve as the natural watershed between the Yangtze River basin and the coastal river basins, and are usually perceived as an important zoogeographic boundary (Zhang and Chen, 1997; Tang W. Q. et al., 2022) and barriers to dispersal for many fish species (e.g., Cao et al., 2013; Chiang et al., 2013; Zheng et al., 2015). Geological evidence have indicated that the large-scale uplifts of the Huangshan Mountains and the Wuyi Mountains took place during the Quaternary (Chao, 1981; Zhang et al., 1990; Wang and Zhu, 1997; Huang et al., 1998). This period correlated with the rapid uplifting stages of the Qinghai-Tibetan Plateau (QTP) since the late Pliocene at around 3.6 Ma (Li and Fang, 1999). The intensive uplifting of the QTP has been assumed to act as a dominant driving force of tectonic movements in eastern China (Xu and Zhou, 2012), resulting in a terrain gradually sloping from west to east and new landscape features (Shi et al., 1998; Zhang et al., 2000; Yan Y. et al., 2018). The early diversification of lineage A, B, and C could be considered as eastward-westward divergence of the ancestral population. The temporal diversification was dated back to the late Pliocene and early Pleistocene (3.00, 2.61, and 2.12 Ma), largely coincide with the orogeny of mountains in southeastern China (like Wuyi Mountains and Huangshan Mountains) and the eastern margin of the QTP during this period. Despite that lineage A is locally endemic to a small river basin in the present, it could be once widespread but only eventually survived in this small region after its isolation. The uplift of mountain ranges likely acted as a natural divide between formerly connected rivers, impeding dispersal and gene flow and resulting in the isolation of populations on both slopes. As a small resident stream fish, the pale chub might be prone to the alternation in landscapes and water systems, resulting in the formation of current phylogeographical structures. Recent phylogeographic studies on two gudgeon species have congruently revealed similar scenarios, in which the significant genetic differences between the coastal river populations and Yangtze River populations were attributed to the long-term isolation between these river drainages as a result of the orogenic events in southeastern China (Yang et al., 2022; Li M. et al., 2023). These findings, along with ours, emphasized the important role of the low mountains and hills in structuring the phylogeographic patterns of small freshwater fishes in southeastern China.

4.1.2. Pleistocene climatic fluctuations facilitated colonization and divergence along the coast

The modern lineage D, E and F occupy the southeastern coastal area of China, and a southward diversification pattern along the coast could be clearly observed. The divergence time between lineage D and lineage E was estimated to be 1.33 Ma and the split between lineage E and F were estimated to be 0.95 Ma, both correspond to the timing of the Early-Middle Pleistocene transition (EMPT, approximately 0.8 ~ 1.2 Ma). In this period an intensification in amplitude and duration of climatic cycles marked the onset of a pronounced succession of glacial and interglacial which characterized the late Pleistocene (e.g., Clark et al., 2006; Ehlers et al., 2018). Geological evidence has suggested that during the EMPT, 25-50m and 20-30m changes in the sea-level took place between 1.2 ~ 2.8 Ma and 0.9 ~ 1.0 Ma, respectively, (Kitamura and Kawagoe, 2006; Berends et al., 2021). The amplitude in eustatic sea level changes of the EMPT were incomparable to those of the LGM, which were 120–130 m lower than the present (Lambeck and Chappell, 2001). However, given that the water depth of inner continental shelf in the East China Sea was usually lower than 60m (Yi et al., 2014), one could still suspect that the alternations in water course due to sea level changes could exert a great impact on the evolution of freshwater fish. The coastal area was impacted by lower sea levels during glacial periods, exposing large areas of continental shelf and paleo-rivers. These paleo-rivers likely increased population connectivity through river anastomoses, which would promote dispersal and gene flow and result in the colonization into new habitats. During the subsequent interglacial periods, sea-level highstands might induce a fragmentation of the formerly connected river systems and result in the isolation of populations. The cyclic changes in water courses could have enabled the southward colonization of the ancestral population and independent evolution after isolations. Drainage alternations associated with sea-level fluctuations during the Quaternary may be of great importance to the range dynamics of freshwater species in coastal China, as many previous phylogeographic studies have suggested (Yang and He, 2008; Yang, 2012; Yu et al., 2014; Yang et al., 2016).

Our findings also shed lights on the origin of primary freshwater fish fauna of the Taiwan Island. Previous studies had widely assumed that rivers in Zhejiang and Fujian region (e.g., Qiantangjiang River and Minjiang River) possibly act as major dispersal routes from continental East Asia to the Taiwan Island, especially for freshwater fishes that restricted to the northern or northwestern part of the island (e.g., Wang et al., 2004; Xu et al., 2005; Lin et al., 2010; Yang, 2012; Ju et al., 2018). The natural distribution of O. acutipinnis-evolans complex in Taiwan was confined to the northern part of the island (Lin, 2008). Our results agree fairly well with this hypothesis. The phylogeographic analysis indicated that populations in Taiwan derived from the Zhejiang-Fujian region via two independent colonization events in different period. Samples from the northern part of the island, including the Danshuihe River and Shuangxi River, showed a close relationship with those sampled from the Qiangtangjiang/ Oujiang River (Figure 1D), while samples from northwestern part of the island derived from those sampled from Mulanxi River (Figure 1F). Geological and phylogeographical evidence have shown that a land bridge recurrently connected the Taiwan Island with the continental East Asia during the Quaternary glacial cycles, facilitating biotic interchanges between both sides (e.g., Wang et al., 2007; Chiang et al., 2010; Lin et al., 2022). The periodic emergence of the land bridge could be a plausible explanation for the multiple colonization routes at different time stages as some previous studies suggested (Yu, 1995; Wang et al., 2004; Chang et al., 2016).

4.2. Genetic diversity

The populations of O. acutipinnis-evolans complex in southeastern China exhibit a high level of genetic diversity. Over 80% of the haplotypes were unique to specific populations, which suggests the possibility of population isolation. This phenomenon may be attributed not only to geographic factors, but also to intrinsic properties of this group of fishes. The O. acutipinnis-evolans complex is a group of small resident fish that prefers lotic habitats, typically found in small montane streams and tributaries. The low dispersal potentials and ecological preference of these fishes could render them vulnerable to the alternation in landscape and hydrological patterns, leading to population isolation or restricted gene flow among populations. Similar patterns have also been observed in other small freshwater fishes in southeastern China (e.g., Yu et al., 2014; Yang et al., 2022).

Grant and Bowen (1998) have proposed four categories of population histories based on haplotype and nucleotide diversity. The first category includes populations with both low haplotype and nucleotide diversity (h < 0.5, π < 0.5%), which may have underwent recent population bottleneck or founder event. The second category comprises population with high haplotype diversity (h > 0.5) but low nucleotide diversity (π < 0.5%), which could be attributed to rapid population growth after population bottleneck. The third category represents populations that have low haplotype diversity (h < 0.5) but high nucleotide diversity (π > 0.5%), which may have experienced bottleneck in a previously large and stable population, or secondary contact. The forth category includes populations that have both large haplotype and nucleotide diversity (h > 0.5, π > 0.5%), which is likely associated with large, stable populations or secondary contact of previously isolated populations. The nucleotide and haplotype diversity were both high (h = 0.9655 ± 0.0028, π = 0.157542 ± 0.076141) in pooled and most individual populations of the O. acutipinnis-evolans complex, confirming the classification in category four. These fishes are typically found in large, stable populations. In general, effective population size and molecular diversity are positively correlated, as observed in many freshwater fishes (e.g., Chiang et al., 2013; Wang et al., 2022).

According to Grant and Bowen’s study (1998), geographic co-occurrence between previously allopatric lineages could lead to large nucleotide diversity. Our study found that several coastal river populations (population 13, FD; population 23, Bs; population 24, ZP; population 25, HA) exhibited higher nucleotide diversity than other populations due to the admixture of two major lineages (Figure 1A), while the remaining populations had a single lineage. Some previous studies focusing on both animals (e.g., Zhao et al., 2008; Li M. et al., 2023) and plants (e.g., Nettel et al., 2008; Havrdová et al., 2015) have attributed the elevated nucleotide diversity to a result of secondary contact. In our case, population 13 (FD) was originated from a small coastal stream (Tongshanxi River). It is plausible to attribute this phenomenon to secondary contact that arise as the consequences of the drainage alternation due to sea level changes or river capture events. Though the relative role of them requires further investigation. Population 23 (BS), population 24 (ZP) and population 25 (HA) all contained haplotypes from both lineage E and F. These populations were all derived from upper or middle reaches of the Jiulongjiang River. Since the lower reach population 26 (CTi) did not exhibit admixture, secondary contact via river capture could possibly contribute to this scenario.

4.3. Deep genetic differentiation and notes on taxonomy

Phylogeographic analysis can provide insights into the speciation processes as it elucidates the presence of incipient subpopulations on the way to becoming new species or that may already represent cryptic species (e.g., Marshall et al., 2011). Our Cyt b gene tree revealed six distinct, almost allopatric evolutionary lineages. Large genetic distances can be observed between some lineages, which could be considered as subspecies or even species level divergence as some previous studies have suggested (e.g., Yu et al., 2014; Chen et al., 2017; Li et al., 2020). However, our morphological analysis concerning both morphometrics and meristic counts failed to distinguish all the six major mitochondrial lineages recovered by the Cyt b gene tree, suggesting that O. acutipinnis-evolans complex could be morphologically conserved. The presence of cryptic species has been suggested to explain the high intraspecific divergence in the opsariichthyine fishes (Wang et al., 2019). Morphologically cryptic lineages could have different adaptations in other aspects, such as behaviors, ecology and physiology (e.g., Braune et al., 2008; Feckler et al., 2014; Taylor and Friesen, 2017). For example, the allochrony in migration and breeding has been assumed to act as an prezygotic isolation mechanism in maintaining the cryptic diversity of the morphologically conserved sand martin subspecies (Tang Q. et al., 2022). Therefore, we still cannot resist that cryptic species might occur. Given the genetic distinctness and the morphological conservatism among lineages, it is still more appropriate to consider this group of fishes as a species complex.

According to Huynh and Chen (2013) and Wang et al. (2019), the two currently described species, O. acutipinnis and O. evolans, can be distinguished by: (1) two distinct mitochondrial lineages with genetic divergence of 2.5%–2.7% in D-loop sequences or 5.2% in Cyt b sequences; (2) allopatric distribution, with O. acutipinnis occupies the Yangtze River basin and Pearl River basin and O. evolans inhabits coastal drainages in southeastern China; and (3) variation of a single branched pectoral fin ray and a single scale above the lateral line, differences in predorsal scales (18–20 in O. acutipinnis vs. 16–17 in O. evolans), and the relative position of the posterior end of maxilla and pectoral fin. Compared to the previous works, our Cyt b gene tree did recovered two monophyletic groups that correspond to these two previously defined species: lineage D + E + F correspond to O. evolans and lineage B corresponds to O. acutipinnis. However, the results of morphological examination disagree with the previous studies that these lineages could be separated by the aforementioned morphological features. The differences of merely a single branched pectoral fin rays, a single scale above the lateral line can be attributed to intra-population variations. They might not be appropriate to act as diagnostic features for distinct species. The length of pectoral fin varied between individuals, age, and sex, with longer pectoral fin that extending way beyond ventral fin origin observed mostly in sexually active males. Therefore, distinguishing O. acutipinnis and O. evolans with the previously proposed morphological diagnostic features would be problematic. Besides, the results of our phylogeographic analysis detected lineage A from the coastal populations and lineage C out of the Yangtze populations. Arbitrarily assigning the coastal lineages as O. evolans and Yangtze-Pearl lineages as O. acutipinnis would also be a hasty decision. Thus distribution as an identification key for both species should also be used with care. Since morphological characteristics and distribution cannot effectively distinguish these two species, DNA barcoding becomes a safer way in their identification. Still, the taxonomy of the O. acutipinnis-evolans complex should be revised carefully, addressing not only morphological differences but also other aspects like ecology and behavior. In addition, the present work, as well as the former ones, only used mitochondrial markers to infer the evolutionary history of the O. acutipinnis-evolans complex. The mitochondrial DNA can only represent a tiny fraction of the whole genome and its limitations have long been addressed (e.g., Galtier et al., 2009). In future studies, it is recommended that additional data, especially from nuclear genome, be incorporated to produce more robust results.

4.4. Implications for conservation and management

Evolutionary significant units (ESUs) and management units (MUs) are important concepts in conservation genetics (e.g., Moritz, 1994; Fraser and Bernatchez, 2001). According to Avise (2005), an evolutionary significant unit (ESU) is one or a set of conspecific populations with relatively long term evolutionary history mostly separate from other such units. Our phylogeographic analysis demonstrated six allopatrically distributed lineages that underwent long geographical isolation. Given this phylogeographic distinctiveness and evolutionary independence, each major lineage should be regarded as a different evolutionary significant unit (ESU). Neglecting the uniqueness of different lineages will certainly increase the chance of genetic homogenization and contamination of the gene pool, which may increase the risk of extinction of local populations. The aforementioned Chinese giant salamander serves as a tragic example (Yan F. et al., 2018). Previous research on the Japanese pale chub Zacco platypus revealed that a distinct lineage native to the Lake Biwa was introduced to other locations in Japan, and resulted in the successful establishment of this lineage in the new environments and widely hybridization with the endemic populations throughout most of the areas in the Japanese islands (Takamura and Nakahara, 2015). Similar situations of other fishes also call for attention to the phylogeographic structures and genetic uniqueness of populations (e.g., Jang-Liaw et al., 2019; Watanabe et al., 2020). Thus, genetic and phylogeographic assessment are crucial for the seemingly widespread but morphologically conserved species, especially when it comes to implementing conservation actions such as translocation or reintroduction.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Ethics statement

The animal study was reviewed and approved by Animal Care and Use Committee of Institute of Hydrobiology, Chinese Academy of Sciences.

Author contributions

HL conceived the research and provided financial support. JG collected samples, conducted laboratory work, and performed data analysis. DY contributed to the sampling and data analysis. JG wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by the National Key R & D Program of China (2018YFD0900806), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB31000000), and National Natural Science Foundation of China (31872234).

Acknowledgments

We thank Siqing Liu, Huaming Hu, Mingzheng Li, and other colleagues for their assistance with specimen collection. We also thank Ping Yang, Lan Zhu, Xianglong Jiang, Qiaoling Deng, Tingyue Yang, Xiaojuan Zhao, and Yahan Wu for their help with laboratory work and data analysis. Thanks also to Tereza Mchura Magati for improving the languages.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2023.1142810/full#supplementary-material

References