fishes Otolith Microchemistry and Demographic History Provide New Insight into the Migratory Behavior and Heterogeneous Genetic Divergence of Coilia grayii in the Pearl River

: Coilia grayii is the anadromous form of anchovy that is distributed in the East and South China Seas. It is a common ﬁsh species in the estuarine area of the Pearl River. Nevertheless, freshwater populations appear upstream in the Pearl River, but the migratory pathway has been mostly impeded by dam construction. Behavioral differences and constrained habitat within tributaries are suspected of promoting genetic divergence in these populations. In this study, we investigated the migratory behavior and genetic divergence of six populations of C. grayii fragmented by dams based on the otolith strontium/calcium (Sr/Ca) ratio, mitochondrial DNA, and microsatellite genotyping. All populations were in freshwater with low Sr/Ca ratios, except the estuarine population (Humen population) hatched in brackish water. Reduced nucleotide diversity corresponding to distance was observed. Populations from distant hydrological regions exhibited a decline in genetic diversity and a signiﬁcant difference with the remaining populations after ﬁtting the isolation by distance model. Pairwise ﬁxation indices conﬁrmed these results and moderate and signiﬁcant differentiation was found between Hengxian site and downstream sites. Furthermore, STRUCTURE analyses revealed that all separated populations exhibited an admixed phylogenetic pattern except for individuals from the Hengxian locality. The upstream sites showed signiﬁcantly increased resistance to gene ﬂow from the estuarine population because of isolation by the dam. The results of the neutrality test and Bayesian skyline plots demonstrated complex demography—individuals’ experienced historical expansion and partial upper-dam populations had recently undergone a colonization, forming a new genetic structure. Accordingly, this study demonstrates differences in the migration pattern and genetic differentiation of C. grayii as a consequence of demographic history and current processes (habitat fragmentation and colonization).


Introduction
Improved understanding of anadromous fish population dynamics requires identifying the ecologically different populations, the genetic structure, and a better knowledge of the life histories of individuals in coastal habitats [1,2]. It is difficult to predict the migratory patterns and population diversity of anadromous and estuarine fishes that are not well-resolved, because of variability in movement and great adaptability to the heterogeneous environments [3]. Several organisms even take residence in opposite environments within diverse life histories, potentially derive into subpopulations [4,5], and display population structures that are dependent on the physical limitations of the habitat and their dispersal capability [6], including many anadromous species of salmon [5,7], some species of sturgeon [8,9], Hilsa shad Tenualosa ilisha [10,11]. These has important ramifications for population connectivity, population persistence, and the potential management of fish stocks.
Long-distance dispersal and clear barriers to gene flow affect population connectivity and induce genetic divergence [12,13]. Furthermore, the colonization and isolation of ancestral species can lead to locally adapted lineages associated with environmental heterogeneity [14][15][16]. Thus, strong migratory connectivity of anadromous fish provides an excellent opportunity to study the interaction between long-distance dispersal and geographic barrier. Identifying behavioral polymorphisms and exploring their cause promotes an understanding of how creatures adapt to the brackish water habitat and provides valuable insights for discerning the evolutionary forces on genetic divergence.
Coilia grayii (Clupeiformes: Engraulidae) is a common anchovy that travels long distances in the estuarine area of the Pearl River and is distributed along the coasts from the East China Sea to the South China Sea as a forage fish [17,18]. During the spawning season (March-June), schools of anchovy stocks migrate between their spawning, feeding, and nursery grounds, while traveling in fresh or brackish water [19,20]. Surprisingly, some freshwater anchovies are found in the upper reaches of the Pearl River, which are hundreds of kilometers away from shore. Coilia nasus, a sister taxon of C. grayii, has various migratory behaviors, including freshwater residence, estuarine originating, anadromous, and non-anadromous [21,22]. A recent study revealed the migratory behavior of C. grayii in the Pearl River estuary based on the otolith strontium/calcium (Sr/Ca) ratio using an electron probe microanalysis (EPMA) and indicated that the fish hatches in freshwater [18]. However, it is unclear whether this dispersal behavior has a genetic basis.
Discerning the migratory behavior and genetic structure of C. grayii in the Pearl River is key to identifying phenotypic and genetic polymorphisms and determining the differentiation of the population. We suspect that C. grayii is subject to different evolutionary forces across space within its wide distribution range. For instance, dispersal and isolation have potential consequences for gene flow and exhibit a central role in life-history evolution. The genetic structure of a species with a wide habitat range may exhibit a pattern of isolation by distance (IBD) and reduced gene flow due to restricted movement within rivers. Compared to the estuarine population, geographically distant populations in freshwater have different life histories, and the process of increasing genetic differentiation is commonly assumed to be correlated with increasing geographic distance. If dispersal has a genetic basis, evolutionary forces on the dispersal process may produce predictable differences in dispersal phenotypes among different populations. Additionally, anthropogenic barriers, such as dams, impede upstream migration and produce the founder effect, if some individuals can pass through them, which may lead to another life-history style [23][24][25]. Several hydroelectric dams have been built in the tributaries of the Pearl River during the few last decades, which impede the migratory pathway of long-distance diadromous fish, reduce gene flow, and increase genetic drift [26][27][28]. Fortunately, over the last decade, rapid technological development in the field of genetics has allowed the study of dispersal pathways and population connectivity of long-distance migrants [29]. Diverse genetic markers, such as mitochondrial and microsatellite genotypes, are some of the most robust approaches to disentangling population dispersal patterns [30,31]. For example, the high mutation rate of the mitochondrial D-loop region allows more rapid accumulation of genetic differentiation between local populations and microsatellites serve as a good complement to mitochondrial markers as they have biparental inheritance, thus their combination resulting in better resolution when attempting to differentiate populations [32,33]. In addition, the high degree of resolution better identifies close relatives from different areas that could presumably reflect dispersal processes with potential dispersal events among populations and population structure within a species [34,35].
C. grayii is suspected to be subject to different evolutionary forces across space, in addition to exhibiting diverse migratory behaviors and life histories that interact to accelerate differentiation. We hypothesized that connectivity was lower within upstream populations of C. grayii than within downstream populations with concomitantly higher levels of genetic differentiation. In this study, we used the otolith Sr/Ca ratio to determine the migratory behaviors of C. grayii in the tributaries and estuary of the Pearl River. A comprehensive multi-locus dataset (cytochrome b gene and the D-loop region of mtDNA, microsatellites fragment length) was utilized to investigate the genetic diversity and population structure. Specifically, we tested whether C. grayii displays strong genetic divergence across a large geographical distance. Finally, we distinguished the genetic structure of the individuals captured from the six sites in three watersheds and determined whether it coincided with the geographic distribution pattern, and inferred the drivers and demographic history.

Sample Collection
A total of 180 specimens were collected from fish docks or local markets in 6 sites of the Pearl River from November 2015 to May 2016, including Qingyuan (QY) and Lixi (LX) in the Bei River tributary; Boluo (BL) and Huizhou (HZ) in the Dong River tributary; Hengxian (HX) in the Yu River tributary and Humen (HM) near the estuary ( Figure 1). They are separated by dam construction at different times (Table 1). We defined Humen (HM) as an estuarine population; populations that are downstream from a dam are downstream populations (BL, QY), the other are upstream populations (HX, LX, HZ). Sagittal otoliths were collected and initially cleaned with 75% ethanol. Standard length and body weight were also measured for each individual. Dorsal muscles were sampled and transferred to the laboratory in liquid nitrogen and then stored at −80 • C for further experiment. different areas that could presumably reflect dispersal processes with potential dispersal events among populations and population structure within a species [34,35]. C. grayii is suspected to be subject to different evolutionary forces across space, in addition to exhibiting diverse migratory behaviors and life histories that interact to accelerate differentiation. We hypothesized that connectivity was lower within upstream populations of C. grayii than within downstream populations with concomitantly higher levels of genetic differentiation. In this study, we used the otolith Sr/Ca ratio to determine the migratory behaviors of C. grayii in the tributaries and estuary of the Pearl River. A comprehensive multi-locus dataset (cytochrome b gene and the D-loop region of mtDNA, microsatellites fragment length) was utilized to investigate the genetic diversity and population structure. Specifically, we tested whether C. grayii displays strong genetic divergence across a large geographical distance. Finally, we distinguished the genetic structure of the individuals captured from the six sites in three watersheds and determined whether it coincided with the geographic distribution pattern, and inferred the drivers and demographic history.

Sample Collection
A total of 180 specimens were collected from fish docks or local markets in 6 sites of the Pearl River from November 2015 to May 2016, including Qingyuan (QY) and Lixi (LX) in the Bei River tributary; Boluo (BL) and Huizhou (HZ) in the Dong River tributary; Hengxian (HX) in the Yu River tributary and Humen (HM) near the estuary ( Figure 1). They are separated by dam construction at different times (Table 1). We defined Humen (HM) as an estuarine population; populations that are downstream from a dam are downstream populations (BL, QY), the other are upstream populations (HX, LX, HZ). Sagittal otoliths were collected and initially cleaned with 75% ethanol. Standard length and body weight were also measured for each individual. Dorsal muscles were sampled and transferred to the laboratory in liquid nitrogen and then stored at −80 °C for further experiment.

Otolith Preparation and Microprobe Analysis
Only left sagittal otoliths of 16 specimens (every 3 individuals from the BL, the LX, the QY, the HM; every 2 individuals from the HX and the QY) were selected for otolith analysis whose body length ranges from 172 to 240 mm. The otoliths were cleaned and air-dried and embedded in epoxy resin (Epofix, Struers). The otoliths were then grounded to expose the core in the sagittal plane using a grinding machine equipped with a diamond cup wheel and polished further with 6 and 1 mm diamond paste on an automated polishing wheel. Finally, the otoliths were cleaned in an ultrasonic bath using distilled water, dried at room temperature, and coated with carbon for EPMA (electron probe microanalysis) measurement.
X-ray intensity maps of strontium (Sr) in the sagittal plane of the otoliths of 16 individuals were performed for "life history transect" analyses of Sr and Ca concentrations along the longest axis from the otolith core to the edge using EPMA (JEOL JXA-8230) as described in Dou and Yokkaichi (2012) and Jiang [22]. Calcite (CaCO 3 ) and strontianite (SrCO 3 ) were used as standards. The accelerating voltage and beam current for X-ray intensity maps were 15 kV and 5.0 × 10 −7 A, respectively, and the electron beam was focused on a point 10 µm, pixel size was 7 × 7 µm and counting time was 50 ms. The beam current changed to 2.0 × 10 −8 A to measure the concentrations of Sr and Ca. The electron beam was focused on point 5 µm with measurements spaced at 15 µm and the counting time for each point was 15 s. Statistical analysis was performed in R. Sr/Ca concentration ratios were traditionally determined and represented as Sr/Ca × 1000. The Sr/Ca ratios and X-ray intensity maps of Sr content were used to define the habitat phases in the otolith by comparing them to known patterns of various salinity exposure phases in C. grayii and C. nasus [22,36,37], where a Sr/Ca ratio of ≤3 is considered to be indicative of a freshwater phase and is colored blue; a ratio of 3-5 is considered indicative of a brackish water phase and is colored green-yellow, and a ratio > 5 is considered indicative of a fully seawater phase and is colored red. The color map patterns of otolith Sr were employed to provide more objective and intuitive insights regarding fish migratory patterns and environmental salinity history.
Amplification reactions of the mtDNA were set up in a total of 25 µL volumes, and amplification reaction mixtures consisted of 100 ng DNA template, 12.5 µL 2× Easy Taq ® PCR SuperMix (Transgene, Beijing, China), 1 µL of 2.5 µM primers each, with sterilized water added to make up the final volume to 25µL. The optimal conditions were: 95 • C for 5 min, 35 cycles of 95 • C for 50 s, 55 • C for 40 s, and 5 min at 72 • C. The amplification prod- ucts were analyzed by 1.0% agarose gel electrophoresis. All positive PCR products were sent to the Beijing Genomics Institute (BGI) for sequencing. At last, we got 155 sequences for D-loop and 132 sequences for cytochrome b (cyt b).
34 microsatellite primer pairs previously isolated from Coilia nasus [39] (Fang et al., 2012), of these, only 5 primer pairs amplified fragments of the expected size, which were genotyped for 174 individuals following the protocol of YangQ et al. [21] (2010) ( Table S1). Amplification reactions of the microsatellites were also set up in a total of 25 µL volumes, consisting of 12.5 µL 12.5 µL 2× Easy Taq ® PCR SuperMix (Transgene, Beijing, China), 1 µL of 5 µM the Tail A (FAM), 0.8 µL of 5 µM the forward primer, 2.5 µL of 5 µM the reverse primer, 25 ng of total DNA and a variable amount of RNase-free water to reach the final volume. All reactions were run in a TGradient thermocycler as follows: 95 for 15 min, 35 cycles of 94 for 30 s, 55 for 90 s, and 72 for 60 s, and the final extension at 72 for 10 min. PCR products were analyzed on 6% denaturing polyacrylamide gels using a Sequi-Gen ® GT sequencing apparatus (Bio-Rad, Hercules, CA, USA) with 97 Sharktooth lanes (Bio-Rad), then visualized by silver staining, carried by BGI. The results were diagnosed and evaluated using the GeneMarker v2.2.0 software.

Data Analysis
The mtDNA sequences were aligned and checked visually in the software Seqman II (DNASTAR) and MEGA 6.0 [40], then refined with Gblock v0.91b [41,42] using default gap penalties. Modeltest v2.3 [43] (Nylander 2004) was used to identify the appropriate model of sequence evolution. Mitochondrial DNA diversity was evaluated using haplotype diversity (Hd) and nucleotide diversity (π) for each population using DnaSP v5.0 [44]. A median-joining network of mtDNA control region and cyt b [45] was estimated in PopART (http://popart.otago.ac.nz, accessed on 12 November 2018). Using the program NETWORK v4.6.1.0 (copyright 2004-2012, Fluxus Technologies Ltd., Colchester, UK; http:// www.fluxus-engineering.com, accessed on 12 November 2018), median-joining haplotype networks were drawn to visually illustrate haplotype variability. Evolutionary relationships among haplotypes were reconstructed by Bayesian inference in MrBayes 3.2 [46] with the substitution model selected by the Modeltest 2.3 and mtDNA sequence of C. nasus (NCBI NO. AP009135) as a close relative outgroup. We implemented MrBayes with two independent parallel runs, four incrementally heated Metropolis-coupled Markov chain Monte Carlo (MCMC) in each run, ran for more than 5 million generations with burn-in 1 million until the average standard deviation of split frequency was below 0.01.
Pairwise genetic distances between populations or individuals were estimated using fixation indices (F ST , Φ ST ). The significance of the fixation index was tested by 1000 permutations in Arlequin v3.5. A hierarchical analysis of molecular variance (AMOVA) was used to estimate the differentiation at the population level and the individual level. In addition, we also tested for isolation by distance (IBD) by implementing the Mantel Test [47], testing the correlation between the geographic distance and genetic distance (the F ST distance matrix based on the 5 microsatellite loci and Φ ST distance matrix based on D-loop sequences). The genetic distance was measured using , while the geographical distance was estimated from corresponding coordinate distances of all sites along the surface of the earth. All the analyses were performed in the R packages "SoDA" and "vegan".
To investigate past demographic expansion, two neutrality tests, Tajima's D [48] and Fu's Fs [49], as well as mismatch distributions [50] were implemented in DnaSP version 5.1 [44]. The Bayesian skyline plot was used to infer the historical population dynamics of Coilia grayii based on the sampled sequences of the control region. It was computed in BEAST v 2.5.1 [51] with a chain length of 100 million generations and a sampling interval of 10,000. The chain convergence and the skyline plot graphic were visualized in the software Tracer v1.5 [52], keeping the effective sample size (ESS) higher than 200. We performed the EBSP-Extended Bayesian Skyline Plot method [53]. Based on microsatellite genotypes, we used Popgene v1.32 and GENALEX v6.5 to assess deviation from Hardy-Weinberg equilibrium (HWE), including polymorphism information content (PIC), observed heterozygosity (H o ), expected heterozygosity (H e ), the heterozygote deficit within populations (F IS ) for each population, and AMOVA. Pair F ST was calculated by Arlequin v3.5, and Nm is inferred from [54] was also used to test population structures based on microsatellites. The iterations of the burn-in period and the number of Monte Carlo Markov chain (MCMC) runs were fixed at 200,000 and 1,000,000 under the admixture model and the assumption of correlated allele frequencies among populations. We conducted five replicated runs for each cluster value of K from 1 to 6. The results from STRUCTURE were uploaded on the website of STRUCTURE HARVESTER based on the ∆K method [55]. After defining K's value, we used CLUMPP v1.  (Piry et al., 1999) was used to assess possible deviations in mutation-drift equilibrium in C. grayii. To analyze the data set of microsatellite genotypes, we used the sign test and Wilcoxon signed-rank test with the infinite allele (IAM), the stepwise mutation (SMM), and two-phase mutation (TPM) models (90% stepwise and 10% multistep mutations 1000 iterations).
The otolith microchemistry of Humen was indicative of the different migration patterns (estuarine migratory). The Sr/Ca ratios of individuals examined from the otolith core to the outermost regions (HM01, HM02, and HM08) were 2.44 ± 1.19, 2.54 ± 1.20, and 3.82 ± 1.53 (Figure 3), which were higher than those collected in freshwater. The X-ray intensity maps revealed that all estuarine otoliths had yellowish or reddish areas in the core region (higher Sr: 3.91 ± 0.69, 4.04 ± 0.85, and 5.11 ± 1.20 for HM01, HM08, and HM02, respectively) as a high signal recorded in Figure 2 and the Sr/Ca ratio subsequently dropped beginning at a range of about 300-500 µm away from the otolith core. These core markers indicated that individuals migrated in chemically distinct environments. The peak Sr/Ca ratio values in the life-history transects of the three estuarine samples were in the core region, which was consistent with migrating between freshwater and seawater, and there were fewer fluctuations of Sr concentration in the peripheral/outermost regions. Marked Sr deposition co-occurred with rapid growth during the early life history of C. grayii with migration into and residency in saltwater (Humen). in the levels of genetic diversity. High haplotype diversity (Hd) in combination with low nucleotide diversity (π) (HdD-loop = 0.853 ± 0.000284, πD-loop = 0.0039 ± 0.00011) was derived from numerous closely associated haplotypes. Among the six populations, the highest Hd and π were recorded in the HM and QY populations, while the HX population presented the lowest Hd and π. For the cytochrome b gene (cyt b), the sequence data were obtained from 132 individuals of C. grayii, which yielded 18 haplotypes and segregated 22 sites (Table S4). The Hd and π values were lower compared to the D-loop (Hdcyt b = 0.306 ± 0.000284, πcyt b = 0.00051 ± 0.00011).

Phylogenetic Analyses
A 1091-bp segment of the mitochondrial control region was amplified from 155 individuals for the six localities. This process produced 40 segregating sites and 37 haplotypes (Table S2). Three dominant haplotypes were shared by four or five populations, while seven haplotypes were distributed in two or three populations, and the remaining haplotypes were single haplotypes (Table S3). Only two haplotypes were detected in the HX population in Yu River and only two dominant haplotypes were found in the LX population of the Bei River upstream. Spatial variations (including long-distance) were observed in the levels of genetic diversity. High haplotype diversity (Hd) in combination with low nucleotide diversity (π) (Hd D-loop = 0.853 ± 0.000284, π D-loop = 0.0039 ± 0.00011) was derived from numerous closely associated haplotypes. Among the six populations, the highest Hd and π were recorded in the HM and QY populations, while the HX population presented the lowest Hd and π. For the cytochrome b gene (cyt b), the sequence data were obtained from 132 individuals of C. grayii, which yielded 18 haplotypes and segregated 22 sites (Table S4). The Hd and π values were lower compared to the D-loop (Hdcyt b = 0.306 ± 0.000284, πcyt b = 0.00051 ± 0.00011).
Based on genetic distances, the Neighbor-Net networks within the D-loop exhibited all of the predicted linking haplotypes in C. grayii, and many less frequently occurring haplotypes differed from the abundant haplotype by a single nucleotide, particularly in the HM assemblages, whereas the BL and QY populations had haplotypes with less frequency and multiple mutations, with few haplotypes among the LX and HX populations (Figure 4b). The cyt b network displayed only 18 haplotypes and one distinct dominant haplotype in the whole population, with several low-frequency haplotypes joining the dominant haplotype (Figure 4a). Based on genetic distances, the Neighbor-Net networks within the D-loop exhibited all of the predicted linking haplotypes in C. grayii, and many less frequently occurring haplotypes differed from the abundant haplotype by a single nucleotide, particularly in the HM assemblages, whereas the BL and QY populations had haplotypes with less frequency and multiple mutations, with few haplotypes among the LX and HX populations (Figure 4b). The cyt b network displayed only 18 haplotypes and one distinct dominant haplotype in the whole population, with several low-frequency haplotypes joining the dominant haplotype (Figure 4a). In addition, five simple sequence repeat markers were amplified in 174 specimens, ranging from 7 to 17 alleles with an average of 13.8 per locus (Table 2). We determined that a large number of homozygotes were apparent in DJ360, which also exhibited the highest F IS value (F IS = 0.4378) and likely suffered allelic dropout. Much lower observed heterozygosity (H O ) was detected in samples from two upstream populations (LX and HX) than other populations. These results indicate the upstream populations were in a departure from Hardy-Weinberg equilibrium (HWE) and the possibility of inbreeding, although many genetic polymorphisms (PIC > 0.5) were observed in all populations.  In addition, five simple sequence repeat markers were amplified in 174 specimens, ranging from 7 to 17 alleles with an average of 13.8 per locus (Table 2). We determined that a large number of homozygotes were apparent in DJ360, which also exhibited the highest FIS value (FIS = 0.4378) and likely suffered allelic dropout. Much lower observed heterozygosity (HO) was detected in samples from two upstream populations (LX and HX) than other populations. These results indicate the upstream populations were in a departure from Hardy-Weinberg equilibrium (HWE) and the possibility of inbreeding, although many genetic polymorphisms (PIC > 0.5) were observed in all populations.

Population Genetics and Structure
Analysis of molecular variance based on two mitochondrial genes was employed to show that most of the genetic variability (Table 3) was due to variance within populations rather than among populations. These results were further corroborated by the microsatellite database. The genetic variance based on the cyt b haplotypes was low and not well resolved, nor was there a clear signature of genetic differentiation among the populations. However, the mtDNA control region analyses resulted in high variance among populations, coupled with a significant difference in ΦST, suggesting that genetic variation

Population Genetics and Structure
Analysis of molecular variance based on two mitochondrial genes was employed to show that most of the genetic variability (Table 3) was due to variance within populations rather than among populations. These results were further corroborated by the microsatellite database. The genetic variance based on the cyt b haplotypes was low and not well resolved, nor was there a clear signature of genetic differentiation among the populations. However, the mtDNA control region analyses resulted in high variance among populations, coupled with a significant difference in Φ ST , suggesting that genetic variation between sampling sites accounted for more than 30% of the total variation. Most genetic differentiation was observed in the HX population: its pairwise genetic distance analysis showed at least a moderate differentiation level (Φ ST > 0.25; Wright 1984) compared with the other populations (Table 4). A similar phenomenon of most genetic differentiation within the farthest population (HX) has been reported in microsatellite data (Table S5). The genetic differentiation between HM and the other two upstream populations was higher than the others. The pairwise Φ STD-Loop results also denoted the fact that genetic variation originated from geographical populations among different tributaries and the estuary (Φ STD-Loop = 0.0658-0.5720, p < 0.01).   Figure 5 provides plots of pairwise genetic distance vs. geographic distance to distinguish the correlations between genetic distance and geographic distance among the populations based on the D-loop sequences and microsatellite data. The Mantel test showed that the divergence of the C. grayii geographic populations in the Pearl River was consistent with the geographical isolation model (IBD model). The genetic differentiation among populations was positively correlated with the geographical distance of the pairwise populations (R 2 = 0.5983, p = 0.0028 **) (Figure 5a) or all six populations, excluding Hengxian (R 2 = 0.4916, p = 0.0167 *) (Figure 5c). Furthermore, a strong correlation was detected with the microsatellite data (R 2 = 0.8878, p = 0.0056 **) (Figure 5b). among populations was positively correlated with the geographical distance of the pairwise populations (R 2 = 0.5983, p = 0.0028 **) (Figure 5a) or all six populations, excluding Hengxian (R 2 = 0.4916, p = 0.0167 *) (Figure 5c). Furthermore, a strong correlation was detected with the microsatellite data (R 2 = 0.8878, p = 0.0056 **) (Figure 5b). Bayesian clustering analyses based on microsatellite loci using the ΔK method suggested that the best K value was 3 ( Figure 6). No distinct genetic structure was observed between up and downstream populations except HX. In three cluster assignments, all individuals from the HX population tended to form a distinct cluster well separated from the rest. Bayesian clustering analyses based on microsatellite loci using the ∆K method suggested that the best K value was 3 ( Figure 6). No distinct genetic structure was observed between up and downstream populations except HX. In three cluster assignments, all individuals from the HX population tended to form a distinct cluster well separated from the rest.

Demography History
Fu's Fs-test and Tajima's D-test based on the mtDNA markers (cyt b and D-loop sequences) yielded negative and non-significant values in the majority of the populations, while the total population results were significant because of the population structure. The significantly negative values manifested an excess of rare haplotypes and low-frequency variants among the populations. We performed an analysis using the BOTTLENECK 1.2 software to test possible deviations in mutation-drift equilibrium. The results were significant (p < 0.050) for the IAM on the Wilcoxon test (Table 5) among three upper-dam populations (HX, LX, and HZ). In part agreement with these findings, the allele frequency distribution pattern of HX was inconsistent with a standard L shape (Figure 7). These results suggested that the populations have recently undergone a colonization and founder effect. Further analysis of the population demographic history is shown by the Bayesian skyline plot of the D-loop sequences-the species has experienced significant expansion in the past 20 thousand years ( Figure 8). Thus, the colonization and recent historical population expansion were the main causes of the demographic changes in the overall population.

Demography History
Fu's Fs-test and Tajima's D-test based on the mtDNA markers (cyt b and D-loop sequences) yielded negative and non-significant values in the majority of the populations, while the total population results were significant because of the population structure. The significantly negative values manifested an excess of rare haplotypes and low-frequency variants among the populations. We performed an analysis using the BOTTLE-NECK 1.2 software to test possible deviations in mutation-drift equilibrium. The results were significant (p < 0.050) for the IAM on the Wilcoxon test (Table 5) among three upperdam populations (HX, LX, and HZ). In part agreement with these findings, the allele frequency distribution pattern of HX was inconsistent with a standard L shape (Figure 7). These results suggested that the populations have recently undergone a colonization and founder effect. Further analysis of the population demographic history is shown by the Bayesian skyline plot of the D-loop sequences-the species has experienced significant expansion in the past 20 thousand years ( Figure 8). Thus, the colonization and recent historical population expansion were the main causes of the demographic changes in the overall population.

Migratory Behavioral Variations
The elemental concentrations in otoliths differed significantly among locations. These concentrations are strongly related to the chemical compositions of the environment in which individuals reside during a particular period of time [59,60]. This study offers new insight into the life history of C. grayii and contributes to determining migratory behavior using otolith microchemistry. We discovered two types of life histories using EPMA-a life history from residency in freshwater, and the migratory life history from the residency of the Humen population in the estuary. This phenomenon reveals that C. grayii is capable of migrating among chemically distinct environments, which is particularly important in poorly studied polymorphic migratory behaviors. The results were consistent with previous findings that both migratory and freshwater resident types are present in C. nasus which is the most closely related sister species of C. grayii in Poyang Lake [61]. The Sr/Ca ratio has been used in many studies as a marker to indicate a shift between freshwater and marine environments in anadromous fish and to study population structure and plasticity in different habitats [62]. Most studies have been conducted in nearshore environments, either within estuaries or on shallow coastal reefs [63,64].
It is generally accepted that C. grayii juveniles in the Pearl River Estuary are migratory fish that complete their early life stages in freshwater [22]. The results of the EPMA analysis of sixteen individuals collected in the estuary and three tributaries appeared in two migratory patterns. All tributary samples were freshwater residents (Sr/Ca ratio < 3), suggesting that they may never enter seawater during their early stages, while estuarine samples presented high Sr/Ca ratios (2)(3)(4)(5) and hatched in brackish water. The microchemical traces in the otoliths were consistent with juveniles drifting with the current into the estuary to grow [65,66].
There were two sites with distinct water chemistry that may have produced distinct signatures in the otolith cores of the fish captured there. The distinct chemical "zones" of high and low Sr/Ca ratios found in the estuarine otoliths suggest that individuals in the Humen population (HM) spent multiple growing seasons in several chemically distinct environments, which was behaviorally distinct from the freshwater fish. The Sr/Ca ratios observed in the present study ( Figure 3) varied widely, but are easily explained as the signature of entry into brackish water from freshwater. Further studies with more extensive sampling will be necessary to clarify these life-history details and determine how common these hypothesized behaviors may be within this widespread species.

Genetic Diversity
The high haplotype diversity (Hd) of the mtDNA D-loop marker in the present study could be attributed to active migration of the anadromous anchovy, particularly in the case of HM (higher derived locus frequency). Populations with a larger effective population size result in a greater number of haplotypes compared to those with smaller effective population sizes [67,68]. The high haplotype diversity values were similar to those reported for widely distributed anadromous fish and semi-migratory fish, such as sardines [69] and Atlantic herring Clupea harengus [10,11] By contrast, low nucleotide diversity implies a high number of closely related haplotypes with very minimal differences between them [70]. Although the numbers of haplotypes identified for all species were large, genetic differences between haplotypes were considerably low. This high haplotype diversity with low nucleotide diversity is similar to many marine taxa inferred using mtDNA sequence data [71,72]. When a new population develops from a small group, novel mutations are accumulated rapidly, accompanied by the increase in population size, but the nucleotide variation does not accumulate in the short term, resulting in a group with high haplotype heterogeneity and low nucleotide diversity [73]. This group usually has one or two widely existing central haplotypes, and the remaining branch haplotypes are derived from the central haplotype (high frequency of latest mutation singletons or doubletons with a widely distributed haplotype giving rise to a star-shaped network) [74]. Overall, the genetic diversity of the C. grayii populations was low (π D-loop < 0.005), while total haplotype diversity was high (Hd D-loop > 0.5) indicating rapid population expansion. It was noteworthy that low observed heterozygosity (H O ) was apparent in two upstream populations from microsatellite markers analysis, which indicates the possibility of the presence of inbreeding or Wahlunds' effect from population substructure [75,76]. Combining both results, it is reasonable to speculate that upstream populations are based on a few founders, which often occur at the edge of an expanding species range exhibiting lower genetic diversity [77], particularly dam separated them in this research.

Population Structure and Habitat Fragmentation
A clear geographic structure was inferred based on the genetic data-the distinctive HX population was separated from the remaining populations and the increasing tendency towards genetic differentiation between LX population and the downstream populations, indicating that the edge population of an expanding species presented genetic differentiation. These patterns raised an intriguing question: why was HX isolated from the remaining populations, while the others were not? Taken together, the complex separation suggests that geographic distances are likely to be related to the population genetics of C. grayii, as distinct levels of genetic differentiation were observed among long-distance sites, for instance, the river pathway between Hengxian and Humen is more than 890 km long. Furthermore, HX is not only just the farthest population, based on the map, its geographical distance to the closest population seems to be much larger than the largest distance among the remaining populations. It is generally accepted that a long geographic distance promotes increased genetic differentiation and decreased genetic diversity [78]. The Mantel test showed that the population genetic distance of C. grayii in the Pearl River was consistent with the geographical isolation model (IBD model). Genetic differentiation among the populations was positively correlated with geographical distance. Considering that the HX population was farthest from the other populations, the result excluding Hengxian was also significant (R 2 = 0.4916, p = 0.0167 *). Similar phylogeographic divergence has been reported in other aquatic organisms (i.e., red swamp crayfish and darters) in the Pearl River basin [79,80]. The IBD pattern also agrees well with the life history of this species.
Barriers that strongly impact wildlife migratory behavior, such as dams, are considered as a driver for affecting genetic connectivity within the river network, resulting in altered genetic structure and favoring genetic divergence, and possibly hindering the recovery of populations [25,81,82]. Genetic differentiation is often identified among populations that are up-or downstream of these types of physical barriers for species including Chinook salmon, brook trout, and chum salmon [83,84]. We demonstrated a reduction of nucleotide diversity presented in the upper-dam populations in the study, which might be due to restricted gene flow and genetic drift. We speculate that the level of population differentiation is positively correlated with the amount of time that a dam has blocked migratory access, although no statistical evidence can be presented. The above-dam population in the Yu River has experienced long-term isolation. The dam was built in 1985 and the river had significantly low genetic diversity (π D-loop = 0.00039, Hd D-loop = 0.071) with only two haplotypes and low gene flow. The upstream populations with long forming time properly have undergone a colonization and founder effect, and the effective population size of the isolated population became smaller. Thus, genetic drift and genetic differentiation among populations increased [85]. Interestingly, the HM and HZ populations had higher divergence than QY, although the distance between Humen and Qingyuan is nearly twice as far as that from Humen to Huizhou.
Strong evidence demonstrates that both geographic distance and habitat fragmentation affect spatial genetic variation; however, we cannot distinguish the interaction between the two components, and the presence of barriers may have a considerable impact on IBD patterns [86]. For instance, the obvious genetic differentiation of population BL and HZ separated by the Jiantan dam should have been observed, but the gene flow still was apparent due to the short distance between the two populations. In this research, the founder effect is a force that could not be neglected and is also capable to contribute to the interaction effect. This is a modest dataset in terms of sample sizes, and we need a broader sampling scheme to eliminate spatial correlation in the loci frequency distributions caused by barriers that are under investigation and test the validity of these hypotheses in the future. Another important issue must be pointed out: By establishing a set of microsatellite markers to use in each analysis for the given species, the number of the markers to employ should be sufficient and highly discriminating. Considering the limitation of the low number of microsatellite markers in this research, the lack of resolution found may be a consequence of that.

Fluctuations in Demography and Migration
The low genetic diversity of the isolated populations (HX) from the distant upper-dam was probably a consequence of decreased heterozygosity (Ho), which led to a decrease in population size and changes in the genetic structure of the population [87]. The microsatellite analysis revealed that the farthest upstream population has undergone colonization, supported by strongly negative values from Fu's Fs indices and Tajima's D neutrality test, which were most likely the result of population declines. These populations are distributed in an isolated watershed with few opportunities to disperse across barriers or long distances. This genetic drift has not taken place in all subpopulations. The downstream populations along the Pearl River Estuary exhibited stable size, where diverse haplotypes and high gene flow precludes the possibility of deep population structuring and impedes speciation [19]. For instance, the high genetic diversity of the specimens from HM, as well as the many haplotypes (singletons or doubletons with a widely distributed haplotype giving rise to a star-shaped network) and high polymorphism are consistent with population stability in this region [20]. The population without restricted movement within the river presented unique genotypes and high levels of diversity, more ancient alleles, and many private D-loops. This pattern is supported by the major role of productive estuarine systems. High local gene flow and a lack of genetic structure also coincide with studies on eastern two-winged flyingfish in the Pacific Ocean [88].
Two distinct life history types (anadromous vs. freshwater resident) in the population reminded us to infer the population dynamic history. Notably, non-significant Tajima's D values in the majority of the populations and Bayesian skyline plot provided evidence that the whole population underwent an episode of recent demographic expansion. These findings confirm that the increase in sea level after the Late Pleistocene resulted in greater availability of habitats when favorable conditions emerged. Similar results have been reported for the expansions of other fish, for example, salmon and silversides [17,77]. New habitats with more ecological resources, such as food, encourage more long-distance dispersal and thus favoring demographic expansion [89,90]. As an anadromous or semi-migratory fish, C. grayii might colonize and some populations take up residence in freshwater, even upstream where there are increases in anthropogenic activities, then, the separated populations are affected by reduced gene flow, ultimately resulting in genetic divergence. Undoubtedly, the heterozygosity of migratory behavior originated from these historical events, and severe environmental fluctuations and early founders could increase the evolutionary load.
In summary, these strongly suggest the impact of isolation and demographic history on the genetic structure and polymorphic migratory behaviors of C. grayii populations, which is considered suitable material for founder effect and divergence research.

Conclusions
The combination of otolith microchemistry and genetic analyses suggested anadromous Coilia grayii were observed in the estuarine population which migrated among chemically distinct environments. Furthermore, we also demonstrated that habitat fragmentation and long geographic distance had a significant effect on the genetic diversity of the upstream populations. The isolation driven by behavior and the dam restricted gene flow and reduced genetic diversity. We offer a robust validation that anthropogenic barriers affect migratory behaviors, genetic diversity, population demographics of wild individuals. Here, we provide a better knowledge of the life histories of C. grayii in coastal and freshwater habitats and their fisheries management should be implemented to assure the connectivity of the habitat and the genetic integrity in these C. grayii stocks, building a foundation and guidelines for wild population protection.