A total of 107 specimens of A. normalis was collected from nine localities on Hainan Island and mainland China (Fig. 1, detailed in Table 1). These populations were sorted into five groups based on the landforms. Fishes were collected from field sites with seines and euthanized with MS-222 (Sigma). Samples were fixed and stored in 100% ethanol. Genomic DNA was extracted from muscle tissue using a Genomic DNA Purification Kit (Gentra Systems, Valencia, CA). The entire cyt b gene and control region fragment were amplified by polymerase chain reaction (PCR) using primers from Xiao et al. (2001) and Zhou et al. (2017), respectively. Each 50 µl PCR mixture contained 5 ng template DNA, 5 µl 10x reaction buffer, 5 µl dNTP mix (10 mM), 5 pmol of each primer and 2 U of Taq polymerase (Promega, Madison, WI, U.S.A.). PCR was performed on an MJ Thermal Cycler with the following program: one cycle of denaturation at 95 °C for 4 min, 30 cycles of denaturation at 94 °C for 45 s, annealing at 48 °C for 1 min 15 s and extension at 72 °C for 1 min 30 s, followed by 72 °C extension for 10 min and sample storage at 4 °C. The purified PCR products were sequenced using an ABI 377 automated sequencer (Applied Biosystems, Foster City, CA, U.S.A.). Chromatograms were checked with the Chromas software (Technelysium), and sequences were manually edited using BioEdit 6.0.7 (Hall 1999). Moreover, to examine the taxonomic status among Pararasbora, Yaoshanicus, Nicholsicypris and Aphyocypris, the sequences of Opsariichthyinae from GenBank were obtained (see Results).

Nucleotide sequences were aligned in Clustal X 1.81 (Thompson et al. 1997). Levels of intrapopulation genetic diversity were estimated by indices of haplotype diversity (h) (Nei and Tajima 1983) and nucleotide diversity (θ_π and θ_ω) (Jukes and Cantor 1969) in DnaSP 4.10.8 (Rozas et al. 2003). Comparing estimates generated by these two indices (θ_π, current genetic diversity, and θ_ω, historical diversity) provides insight into population dynamics over recent evolutionary history (Templeton 1993). The existence of phylogeographic structure was examined following Pons and Petit (1996) by calculating two genetic differentiation indices, G_ST and N_ST, in DnaSP. Pairwise F_ST values implemented by DnaSP were used to examine the spatial partitioning of genetic variation among populations. AMOVA (analysis of molecular variance) partitioned were used to the observed variation among samples into within-population (F_ST), within-group (F_SC) and among-group (F_CT) components in Arlequin version 3.5 (Excoffier and Lischer 2010). Arlequin was also used to construct the haplotype network. The program SAMOVA (Dupanloup et al. 2002) was used to identify groups of adjacent sampling locations with the maximum extent of genetic differentiation. These analyses used 500 simulated annealing steps and compared maximum indicators of differentiation (F_CT) when the program was instructed to identify K = 2 to K = 8 partitions of the sampling area.

A phylogenetic tree was created by BEAST 1.8.0 (Drummond et al. 2013), which is a Bayesian statistical framework. Phylogenetic relationships were also inferred using maximum-likelihood (ML) in MEGA 6 (Tamura et al. 2013). The GTR+I+G substitution model was selected using the Akaike information criteria (AICc) in jModelTest 2.0 (Darriba et al. 2012). The time to the most recent common ancestor (T_MRCA) was also calculated with the software package BEAST. A molecular clock was calibrated using a divergence rate of 2% per million years (Yang et al. 2016). We used a strict clock with a Bayesian skyline tree and estimated the divergence times of the major lineages to the most recent common ancestor (T_MRCA). We ran 10⁶ generations. The output was visualized in Tracer v1.6 (Rambaut et al. 2013) to verify that convergence and suitable effective sample size were achieved for all parameters. Burn-in and plots for each analysis were visualized using Tracer. The TREEANNOTATOR in the BEAST package was used to summarize tree data, and the tree was viewed using FigTree v1.3 (Rambaut 2014). For ML analysis, bootstrapping was performed with 1000 replications.

In addition, to determine the possible diversification scenarios of A. normalis, a statistical dispersal-vicariance analysis (S-DIVA), a program that complements DIVA, was employed to determine statistical support for ancestral range reconstructions (Yu et al. 2010). The tree file formats were generated by the program BEAST with 10⁷ Markov Chain Monte Carlo (MCMC) steps and the first 10% as burn-in. The analysis was performed using the ‘maxareas = 5’ option because the populations were divided into six groups (Table 1).

Although the genus Aphyocypris was included in Danioninae, Liao et al. (2011) reassigned this genus into Opsariichthyinae. Thus, sequences of Opsariichthyinae in GenBank from Chiang et al. (2011), Wang et al. (2011) and Lin et al. (2016) were downloaded (Fig. 2). The ML tree (Fig. 2) showed that all sequences fell into four monophyletic groups corresponding to four genera, Aphyocypris, Candidia, Nipponocypris, and Opsariichthys. Among these four genera, the mean pairwise p-distance was 15.84%, with a range of 12.39% to 18.22%. The range of the divergences within these genera was 3.15% to 11.40% (Table 2). In addition, three Aphyocypris species, A. normalis, A. moltrechti, and A. chinensis, are also monophyletic groups. The mean intraspecific divergence within these Aphyocypris species was 1.91%. The mean interspecific p-distance was 9.89%, with a range of 2.63% to 12.59% (Table 2).

Figure 2. 

Phylogenetic relationships of the genera Aphyocypris, Candidia, Nipponocypris, and Opsariichthys using ML analyses of cyt b gene sequence data. Numbers along the branches indicate the percentage of bootstrap support obtained in the ML analyses.

The complete 1,141 base pairs (bp) of cyt b (MH909846-MH909955) and 945 bp of the control region (MH909956-MH910020) sequences were analyzed. A total of 105 haplotypes (2,086 bp, 270 variable sites and 170 phylogenetic informative sites) were obtained for 107 sequences from nine populations analyzed (Table 1; Fig. 1). Nucleotide sequences were A+T rich (64.9%). Haplotype diversities within populations ranged from 0.99 to 1.00 (Table 1). The nucleotide diversity (θ_π) in the region was the highest in northern Hainan Island (average = 1.20) and the lowest in the northern Yunkai Mountains (population HY, 0.23). At the population level, the nucleotide diversity (θ_π) was the highest in population QH (northern Hainan Island; 1.43) and the lowest in population HY (north of Yunkai Mountains; 0.23). Estimates of the current (θ_π) and historical (θ_ω) genetic diversity per site for each sample indicated that all samples showed a pattern of decline (θ_π < θ_ω; Table 1).

All mtDNA haplotypes were population-private haplotypes. Our study found that the sequences were unique, excluding two identical sequences within populations BT and HY. Geographical division assessed by DnaSP indicated significant differentiation among populations (F_ST = 0.530). The pairwise F_ST values ranged from 0.144 (between populations HO and QH) to 0.820 (between populations HY and BT) (Table 3). After pooling by region (Table 1), the pairwise F_ST within the mainland region (mean F_ST = 0.493) was smaller than that within island regions (mean F_ST = 0.514) (Table 3). The mean F_ST between the mainland and island was similar to that within the island (Table 3). A comparison of the fixation indices N_ST (0.532) and G_ST (0.004) revealed that N_ST was much larger than G_ST. This result suggested a strong relationship between phylogeny and geography. These results indicated that the population differentiations were significant.

The Bayesian phylogenetic tree is shown in Figure 3 and is divided into two major phylogroups (I and II). Phylogroup I included seven populations in northern Hainan Island (including) and phylogroup II included only two populations in southern Hainan Island. The best SAMOVA grouping schemes partitioned the sampling area into five groups (K = 5; F_CT = 0.35, p < 0.001). All samples were divided into five units: (1) BS, (2) BT, (3) TH, (4) HY and (5) others. The results of AMOVA indicated significant genetic structures at several levels (Table 4), but the most genetic variability was accounted for by between-group differences among five SAMOVA units. Among these five groups, 34.53%, 24.26% and 41.21% variations existed among groups, among populations within groups and within populations, respectively.

Figure 3. 

BEAST-derived chronograms of 107 mitochondrial DNA sequences of Aphyocypris normalis. The S-DIVA analysis graphical representation of the ancestral distribution is given in the box above the node.

The coalescence time of A. normalis was estimated to be in the early Pleistocene (T_MRCA= 2.33 mya, 1.92 – 2.74). The results of the S-DIVA analysis indicated that possible ancestral populations of A. normalis were distributed on Hainan Island and east of the Leizhou Peninsula. The rising of the WY Range on Hainan Island separated the populations into two groups, south and north of the WY Range. After a vicariance event, the populations in the north of the WY Range migrated to the west of the Leizhou Peninsula.

Today, the genera Pararasbora, Yaoshanicus, and Nicholsicypris are synonymous with the genus Aphyocypris. Our study considered that a genus should fulfill two criteria at least, monophyly and distinctness. In this study, the range of the p-distance among the genera Candidia, Nipponocypris and Opsariichthys was 12.39% to 16.17% (mean = 15.16%) (Table 2). The range of the p-distance between the genus Aphyocypris and the other three genera ranged from 15.31% to 18.22% (mean = 16.52%) (Table 2). Moreover, the mean interspecific p-distance within the genus Aphyocypris was 9.89%, with a range of 2.63% to 12.59% (Table 2). The mean interspecific p-distance within the genus Opsariichthys was 10.69%, with a range of 4.30% to 14.21%. Thus, our study determined that the genera boundary was 15.00%. Based on the genetic variations, our study supported the synonymization of Pararasbora (P. moltrechti), Yaoshanicus (Y. arcus) and Nicholsicypris (N. normalis) with Aphyocypris (Liao et al. 2011). Moreover, Tzeng (1986) indicated that A. normalis is conspecific with A. moltrechti, but our results did not support that finding (Fig. 2; Table 2).

The geological studies proposed that approximately 2–2.5 mya, Hainan was first isolated from the mainland by the results of volcanism and the ground fall in the current Qiongzhou Strait (e.g., Zeng and Zeng 1989). During the Pleistocene glaciations, migrants probably moved between mainland China and the island by the sinking of the Gulf of Tonkin and Qiongzhou Strait (Morton and Blackmore 2001; Yap 2002). Previous studies (Chen et al. 2007; Lin et al. 2016; Yang et al. 2016) have indicated that during the Pleistocene glaciations, the sea level dropped, and the entire region of the Gulf of Tonkin and the Qiongzhou Strait became part of the coastal plain of the Asian continent (Fig. 1). The strait and exposed continental shelves assisted in population dispersions. Thus, the gene flow among the ancestral populations of Hainan Island and the adjacent areas were unlimited. The results presented here suggest that the ancestral populations of A. normalis were distributed widely south of the Yunkai Mountains, including Hainan Island, and east of the Leizhou Peninsula in the early Pleistocene.

Subsequently, our results found that the populations on southern Hainan Island (TH and BT) diverged by a vicariance event. At 0.458 mya, the WY Range arose on Hainan Island and isolated the population in the southern Hainan Island region as lineage II (Fig. 3). This result agrees with previous studies (Yang et al. 2016; Zhou et al. 2017). The WY Range is located in the central and southern regions of Hainan Island and approaches an elevation of 1800 m. The landforms reflect that the rivers on island originate from the central mountainous area and flow outwards. Due to the steep topology of the island, most rivers run directly into the oceans with no connection to the neighboring drainage systems. As in previous phylogeographic studies of freshwater fishes (Yang et al. 2016; Zhou et al. 2017), the WY Range was an important barrier. Zhou et al. (2017) found that the populations of the cyprinid fish O. lepturum on Hainan Island could be divided into three units based on their location in the Qiongzhou Strait, the Gulf of Tonkin and the South China Sea. The genetic structure of A. normalis on Hainan Island showed the same pattern as that of O. lepturum (Zhou et al. 2017). Thus, the WY Range is a phylogeographic break in the phylogeography of A. normalis.

Based on previous studies (Chen et al. 2007; Lin et al. 2016; Yang et al. 2016), the populations on Hainan Island and the mainland diverged by a vicariance event, the Qiongzhou Strait. Previous studies (Chen et al. 2007; Lin et al. 2016; Yang et al. 2016) proposed that populations of freshwater fishes migrated from one side of the exposed Qiongzhou Strait to the other, and the populations of the island and the mainland separated into two highly divergent clades due to the rising sea level. Our study found that the populations of A. normalis on both sides of the Qiongzhou Strait were a mixture (Fig. 3). Although previous studies (Yap 2002; Zhou et al. 2017) considered that the Qiongzhou Strait is a phylogeographic break of freshwater fishes, we found that the populations of G. orientalis (Yang et al. 2016) and O. hainanensis (Lin et al. 2016) on both sides of the Qiongzhou Strait were mixed to some extent. In the population structure of G. orientalis (Yang et al. 2016), lineage I was only distributed to the north of Qiongzhou Strait, but lineage II was distributed both to the north and south of the strait. Likewise, the population structure of O. hainanensis (Lin et al. 2016) displayed the same pattern as that of G. orientalis (Yang et al. 2016). Moreover, geological studies (Voris 2000; Shi et al. 2006; Zhao et al. 2007) suggested that land bridges repeatedly connected the island to the Asian continent. Yang et al. (2016) proposed that when glaciation occurred, the lineages mixed. Thus, the phylogenetic analyses of A. normalis (Fig. 3) in the current study revealed a mixed genetic structure across the strait due to cyclic climate changes. When glaciation occurred again after the WY Range rose, some populations retreated from the mainland to the island. Thus, the phylogenetic analysis (Fig. 3) displayed two divergent phylogroups before the divergences among haplotypes included in the ancestral populations.

Our study found that A. normalis (this study) and G. orientalis (Yang et al. 2016) have similar distribution patterns and that these two species displayed similar phylogeographic patterns. Thus, our study suggested that the effects of environmental changes on Hainan Island and its adjacent area were general patterns. However, the populations of G. orientalis showed moderate to high genetic differentiation, but those within A. normalis showed a high level of genetic differentiation (Table 2). Although these two species both displayed differentiation, there was no shared haplotype between mainland and island populations, and there was no shared haplotype within A. normalis. The haplotype and nucleotide diversities of A. normalis are higher than those of G. orientalis. Moreover, the genus Aphyocypris only includes ten species in East Asia, but the genus Garra contains 138 species in Asia and Africa (Froese and Pauly 2018). Our study found that A. normalis and G. orientalis had different colonization times, with A. normalis in the early Pleistocene and G. orientalis in the late Pleistocene. Our study considered that the ancestral populations of these two species originated from different mainland populations. Thus, our future aim is to understand the comparative phylogeography and geological history in Hainan Island and its adjacent area.