Population structure of Aphyocypris normalis: phylogeography and systematics

Abstract Aphyocypris normalis (Teleostei: Cyprinidae) is an endemic species in South China, but little is known about its genetic structure. This study examined the population structure of A. normalis using sequences of the mitochondrial DNA control region and cytochrome b gene (2,086 bp). In total, 107 specimens were collected from nine populations. All 105 mtDNA haplotypes were identified as belonging to two allopatric phylogroups. The results of a statistical dispersal-vicariance analysis (S-DIVA) suggested that the ancestral populations of A. normalis were distributed widely on Hainan Island and east of the Leizhou Peninsula. A comparison of the fixation indices NST (0.532) and GST (0.004) revealed that the phylogeny and geography had a significant relationship. Our study found that (1) the Wuzhishan and Yinggeling Mountain Range was an important barrier limiting gene exchange between populations on both sides; (2) cyclic climate changes may have shaped migrations and population differentiations; and (3) different colonization times caused different population diversities between codistributed species. In addition, the inter- and intraspecific diversities of the genus Aphyocypris were estimated.


Introduction
Hainan Island is located in the transitional zone between tropical and temperate zones. The island is separated from mainland China to the north by the Qiongzhou Strait and from mainland Vietnam to the west by the Gulf of Tonkin. Hainan Island was first isolated from the mainland by the results of volcanism and rising sea levels approximately 2-2.5 million years ago (mya) (e.g., Zeng and Zeng 1989;Zhao et al. 2007). Based on the ichthyofaunal similarities (Li 1981), Hainan Island was defined as a unit subregion in the South China region. Our study found that some freshwater fishes, e.g., Glyptothorax hainanensis and Opsariichthys hainanensis, are considered endemic to Hainan Island, but previous studies (Chen et al. 2007;Lin et al. 2016) have found that these two species are also distributed in the Pearl River subregion. Moreover, geological evidence indicates that land bridges repeatedly connected the island to the Asian continent during the Pleistocene glaciations (Voris 2000;Zhao et al. 2007). According to phylogeographic studies of freshwater fishes (e.g., genus Glyptothorax see Chen et al. 2007; Garra orientalis see Yang et al. 2016; genus Opsariichthys see , during Pleistocene glaciations, migrants probably moved between mainland China and island populations. The gene flows were interrupted when Hainan Island was separated from mainland China by the sinking of the Qiongzhou Strait (Morton and Blackmore 2001;Yap 2002).
The Wuzhishan and Yinggeling Mountain Range (WY Range) rises steeply from the central and southern regions of Hainan Island and gives way to a broad plain in the north. Lin et al. (2010) found that the WY Range was an important barrier based on the population structure of the Reeves's butterfly lizard. However, Huang et al. (2013) suggested that WY Range did not act as a barrier to gene flow among populations of oriental garden lizards (Calotes versicolor). In addition, the four largest rivers on the island, the Nandu, Changhua, Wanquan, and Linshui Rivers, originate from the WY Range and flow outwards into the Qiongzhou Strait, South China Sea, and Gulf of Tonkin, respectively ( Fig. 1). According to the landforms, the WY Range isolated these four rivers into two regions, southern Hainan (Wanquan and Linshui rivers) and northern Hainan (Nandu and Changhua Rivers). Yang et al. (2016) and Zhou et al. (2017) proposed that the WY Range was an important phylogeographic break based on the population structures of the cyprinid fishes G. orientalis and Onychostoma lepturum. However, Lin et al. (2016) suggested that the WY Range did not act as a barrier to gene flow among populations of another cyprinid, O. hainanensis.
Aphyocypris normalis is a small cyprinid fish restricted to Hainan Island and its adjacent area. Thus, this species is an ideal freshwater fish species to study the biological consequences of the geological history of river systems in this area. Aphyocypris normalis was named Nicholsicypris normalis until 2011. Liao et al. (2011) suggested that Pararasbora, Nicholsicypris, and Yaoshanicus were synonyms of Aphyocypris based on previous morphological studies (Chu 1935;Tzeng 1986;Kottelat 2001). Although Tang et al. (2010) supported their close interrelationship with molecular data, the genetic distances among them were not estimated. Thus, our study examined the taxonomic status of Pararasbora, Nicholsicypris, Yaoshanicus, and Aphyocypris based on molecular data. In addition, the species diversity within Aphyocypris is very low. Liao et al. (2011) found that there are two species within Aphyocypris (A. chinensis and A. kikuchii), and Pararasbra (P. moltrechti), Nicholsicypris (N. normalis), and Yaoshanicus (Y. arcus) are monotypic. Until now, there were eight species within the genus Aphyocypris (Froese and Pauly 2018), and new species were described (Zhu et al. 2013). Thus, our study examined the diversity within Aphyocypris.
The major questions in our study are as follows: (1) What are the taxonomic boundaries within Aphyocypris and its close relatives? (2) How did A. normalis colonize the rivers on the island and mainland? (3) Is there a phylogeographic break in freshwater fish on Hainan Island? To address the aforementioned questions, the mitochondrial DNA cytochrome b gene and control region sequences (hereafter mtDNA) were used to evaluate the phylogenetic relationships and population genetic structure (e.g., Chen et al. 2007;Yang et al. 2012Yang et al. , 2016Chiang et al. 2010Chiang et al. , 2013Hsu et al. 2014). These results will help to establish mtDNA sequence data and develop conservation strategies for freshwater fishes on Hainan Island.

Sampling and molecular methods
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).

Sequence alignment and phylogenetic inference
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 6 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 7 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).

Taxonomic status
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).

Genetic diversity of A. normalis
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).

Population structure and history
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 differen-    = 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.
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.

Inter-and interspecific diversities within Aphyocypris
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).

Phylogeography of Aphyocypris normalis
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  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 ) 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.