The Mitochondrial Markers Provide New Insights Into the Population Demographic History of Coilia nasus With Two Ecotypes (Anadromous and Freshwater)

The Pleistocene climate oscillations have been proven to have profound influences in the evolutionary processes of marine species and left significantly tractable footprints in genomes of marine fauna. However, such history remains unknown for the Japanese grenadier anchovy Coilia nasus, a relatively large fish found in marine, freshwater and brackish water. In order to provide insight into the phylogeographic pattern and historical demography of C. nasus, we employed the mtDNA D-loop and Cyt b as molecular markers to estimate the genetic variation among populations of C. nasus. We took samples from 18 populations throughout the distribution of C. nasus. Partial sequences of mtDNA D-loop (n = 424) and Cyt b (n = 129) were obtained. Four phylogeographic lineages were revealed by the mtDNA D-loop analysis. Frequency changes of the four lineages throughout the species’ distribution indicated limited genetic exchange among three major geographic regions. Patterns of isolation-by-distance were detected among the populations. Neutrality test, mismatch distribution, and network analyses consistently showed that all the four lineages had populations with relatively stable over time and then subjected to sudden expansion during the late-Pleistocene era (c. 20,000–75,000 years ago). The phylogeographic and demographic history of C. nasus was highly influenced by climatic fluctuations in the Pleistocene. Geographic isolation resulted from Pleistocene climatic fluctuations might have profound impacts on the distribution pattern of genetic variation of C. nasus. Our result also indicated that the population expansion event for C. nasus occured along with habitat continuous expansion during the inter-glacial period.


INTRODUCTION
In the Northern Hemisphere, the Pleistocene climate oscillations have been proved to have profound influences in the evolutionary processes of marine species and left significantly tractable footprints in genomes of marine fauna (Hewitt, 1996(Hewitt, , 1999(Hewitt, , 2004Bernatchez and Wilson, 1998;Kitamura et al., 2000;Schmitt, 2007). Representative researches of the origins and evolutions corresponding to glacial refugia and recolonization route have been well studied in many terrestrial and aquatic species (Teacher et al., 2011).
The marine species were always thought to be less sensitive to the effects of Pleistocene climatic oscillations (Hewitt, 2000), however, recent studies have indicated that the Pleistocene glaciations had profound influences on the evolutionary process of many coastal marine species (Avise, 2000;Liu et al., 2007, Wilson andVeraguth, 2010;Liu et al., 2011a). The northwestern (NW) Pacific is characterized by marginal seas and shallow continental shelves, which were deeply impacted by Pleistocene glacial cycles. As one of the largest marginal seas in the northwest Pacific Ocean, most areas of the East China Sea Shelf were exposed during the period of the Pleistocene glaciation maximum (Xu and Oda, 1999). During the last interglacial period, the transgression has happened in the East China Sea Shelf, which caused the coastline to migrate from the Okinawa Trough to the coast of Bohai bay with about 1200 km (Wang, 1999;Xu and Oda, 1999). The sea level rose together to create a large area of newly formed marine habitats, allowing coastal marine species to expand their habitat ranges (Dodson et al., 2007;Liu et al., 2007;Xu et al., 2009;Wilson and Veraguth, 2010). Suitable DNA markers can be used to examine the genetic effects originated from the climatic oscillations. Compared with marine species, anadromous species were more vulnerable to the periodic climatic oscillations in the Northern Hemisphere (Consuegra et al., 2002;Mäkinen and Merilä, 2008;Aldenhoven et al., 2010). For example, recent study shows that the divergence between the resident freshwater and anadromous pairs of the threespine stickleback (Gasterosteus aculeatus) dates to the Late Pleistocene or Early Holocene periods (Bell and Foster, 1994). Phylogeographic studies of the anadromous Atlantic salmon showed that the Baltic Sea basin populations were possibly originated from south-eastern glacial-refugia existed during the Last Glacial Maximum (LGM) period (Tonteri et al., 2005). Historically, the Yangtze River was greatly influenced by the Late Pleistocene climatic oscillations and was likely to be an inland river due to the low sea level, headward erosion and dry climates during the LGM (Yang, 1986;Xiao et al., 2003). The Huai River is much younger than the Yangtze River. The existence of special sediment proves that the Huai River was formed later than 13 ka (Yang et al., 2010a(Yang et al., , 2012. However, a series of geological activities have resulted in the decline of the water level of the Huai River by 13 meters since 2 ka BP (Yang et al., 2010a(Yang et al., , 2012. For these reasons, studies on the phylogeographic and demographic histories of the anadromous species that inhabit Yangtze River, Huai River and the NW Pacific are helpful to illuminate the effects of Pleistocene climatic fluctuations and geographical environment. The Japanese grenadier anchovy, Coilia nasus, is one of the most commercially important fishery resources along the NW Pacific Ocean coasts. C. nasus is widely distributed along the coastal waters of China, the western coast of the Korea, and the Ariake Sea of Japan. The scientific name of this species is always written as C. ectenes. However, as a priority, C. ectenes should be synonymized with C. nasus (Jordan and Seale, 1905). The resources for C. nasus have dramatically declined in the past few years, due to the overexploitation and deterioration of the habitat environment (Cai et al., 1980;Duan et al., 2012).
Coilia nasus has two ecotypes with anadromous form and freshwater form (Cheng, 2011). The anadromous grenadier anchovy spends most of its life in the inshore sea and ascends estuaries and rivers for breeding during its extensive migration, whereas freshwater grenadier anchovy always inhabits freshwater (Yuan et al., 1980;Cheng, 2011). In Japan, C. nasus was only found in the Ariake Sea, which spawned in rivers flowing into the Ariake Bay (Takita, 1967;Whitehead et al., 1988). In China, the anadromous form of C. nasus even penetrated to the Yangtze River for more than 1000 km during the spawning season, whereas the freshwater form inhabited the streams and lakes in the middle reaches of the Yangtze River and the Huai River (Yuan et al., 1976(Yuan et al., , 1980Cheng, 2011). Compared to anadromous form, freshwater form was characterized by earlier age at maturity, shorter maxilla, smaller egg size, lower level of growth rate and smaller body size (Yuan et al., 1976;Ni et al., 1990). Although the anadromous and freshwater forms used to be regarded as different species according to their morphological and ecological differences (Yuan et al., 1976;Cheng, 2011), recent studies demonstrated that the genetic differentiation between the two forms was below the species level (Tang et al., 2007;Yang et al., 2010b;Xu et al., 2020). Until now, enough information about the phylogeographic pattern and historical demography among the freshwater and anadromous populations is missing in the Yangtze River. Previous study on C. nasus has shown that the genetic differentiation between Japanese and Chinese populations may result from life-history traits, ecology factors, and genetic isolation following the glacial retreat (Yang et al., 2011). A deep genetic break was observed among populations of C. nasus based on mtDNA Cyt b information (Ma et al., 2010). The existence of some highly diverged sequences (29 fixed substitutions) in the Fuzhou and Wuhan populations was the main reason for the deep genetic break (Ma et al., 2010).
Population demography and distribution of C. nasus might have been greatly impacted by the Pleistocene climatic oscillations. The anadromous populations of C. nasus prefer estuarine habitats. Most areas of this type are located in the East China Shelf, which was exposed by the lower sea level during the LGM. The Yangtze River, which was also profoundly affected by the lowed sea level in the LGM, is the species' prime spawning area. Thus, we speculate that the phylogeographic and demographic history of C. nasus might be greatly influenced by the Pleistocene climatic fluctuations. Furthermore, the freshwater forms inhabiting in the middle reaches of the Yangtze River might also have originated from the freshwater refugium in the Yangtze River during the last glaciation. We also presumed that the habitat range expansion following the last interglacial period might result in population expansion for C. nasus.
Intraspecific phylogeographical studies have been proven to be important in testing biogeographical hypotheses and making insight into historical events (Bermingham and Martin, 1998).
In addition, few studies have elaborated the genetic structure, demography, and evolutionary process for the freshwater and anadromous populations of C. nasus. The mtDNA is useful for phylogeography and evolutionary studies of animal populations. The mtDNA was characterized with non-recombining rapidly evolving rate, which can be used to detect the past isolations (e.g., in glacial refugia) and secondary contacts from isolated populations (Avise, 2000;Hewitt, 2000). Furthermore, variation in nucleotide level was effective on the inferences of historical demography analysis Ho and Shapiro, 2011). In this study, we analyzed the partial mtDNA D-loop and Cyt b sequences to reveal the phylogeographic pattern and historical demography of C. nasus across its range. The objective of this study was to obtain a comprehensive description of the phylogeographic pattern and historical demography, and then clarify the influence of Pleistocene climate oscillations on the evolutionary history of this species. Also, we aim to determine the glacial refugia for the freshwater populations during the Pleistocene period.

Ethics Statement
All methods were carried out in line with the relevant principle and directions in China.

Sampling and Sequencing
A total of 424 individuals were collected from 18 localities throughout the distribution of C. nasus from 2006 to 2013 (Table 1 and Figure 1). Five specimens of C. grayii were also collected from the Pearl River as the outgroup. All individuals were identified based on morphological features, and a piece of muscle tissue was obtained from each individual and preserved in 95% ethanol or frozen for DNA extraction. Genomic DNA was isolated from the muscle tissue by proteinase K digestion and purified by a standard phenol-chloroform method. The first hypervariable segment of the mtDNA D-loop was amplified with fish primers CR-F (5 -CCACCACTTAGCTCCCAA-3 ) and CR-R (5 -ACCAGATGCCAGGAATAG-3 ). The 5 end of the mtDNA Cyt b gene was amplified for 129 individuals with the primers L14734 (5 -AAC CAC CGT TGT TAT TCA ACT-3 ) and H15149 (5 -CTC AGA ATG ACA TTT GTC CTC A-3 ; Zhang et al., 2007).
The PCR system with 1.25 U Taq DNA polymerase (Takara Co.), 10-100 ng template DNA, 200 nM forward and reverse primers, 200 µM of each dNTP, 10 mM Tris-HCl pH 8.3, 50 mM KCl, 1.5 mM MgCl 2 , was used in the present study. The PCR profile included denaturation at 95 • C for 3 min; 40 cycles of denaturation at 94 • C for 45 s, annealing at 50 • C for45 s, and extension at 72 • C for 45 s; and extension at 72 • C for 5 min in a Biometra thermal cycler. PCR products were purified using the AxyPrep PCR Clean-up Kit [Axygen Biotechnology (Hangzhou) Limited] and were further sequenced on an ABI 3730XL automated sequencer (Sunny Life Technologies Corp., Shanghai, China and Life Technologies, Shanghai, China). Additionally, mtDNA Cyt b sequences for 12 individuals of C. nasus and two individuals of C. grayii were also retrieved from GenBank ( Table 2).

Data Analysis
We used the DNASTAR software (DNASTAR, Inc.) to edit and align all sequences in the study. In order to prevent bias in the estimation of genetic parameters for mtDNA Cyt b, we deleted the tRNA Glu sequence from the 5 end of mtDNA Cyt b fragment. We used the program ARLEQUIN to calculate the population genetic indices such as haplotype diversity (h), nucleotide diversity (π), polymorphic sites, haplotypes distribution, transversions, transitions, and indels (Excoffier and Laval, 2007). We used the software of MrBayes 3.1.2 to construct the phylogenetic relationships among haplotypes (Ronquist and Huelsenbeck, 2003) [44]. We also used the software of MrModeltest 2.3 combined with PAUP 4.0b10 to obtain the best fit model (HKY + I + G) for the MrBayes setting (BIC; Swofford, 2002;Nylander, 2004). Bayesian searches included four chains.
Each chain was run for thirty million generations with a tree sampling frequency of every 100 generations, and the first 25% of the resulting trees were discarded as burn-ins. The haplotype tree was visualized by FigTree (ver.1.3.1; http://tree.bio.ed.ac. uk/software/figtree/). Genealogical structures among haplotypes for the two molecular markers were further evaluated based on the minimum spanning tree (MST) analysis calculated by ARLEQUIN. We also used homologous mtDNA Cyt b sequences of C. nasus and C. grayii to construct the MST to verify the relationships between them (Ma et al., 2010).
We employed the quation µ = dA/2T to estimate the nucleotide substitution rate for grenadier anchovy, where T is an estimate of 3.6 Myr divergence between C. nasus and C. ectenes [this sample should actually be a specimen of C. grayii from Lavoué et al. (2013)], and dA was the net average genetic distance estimated using the homologous sequences of mtDNA D-loop. The genetic distance between lineages originated from the equation dA = dXY -(dX + dY)/2 was calculated based on the Tamura and Nei model using MEGA 3.1 (Kumar et al., 2004). Using this method, the lineage-specific divergence rate for mtDNA D-loop was estimated at about 1.34%/Myr. A 10-fold faster mutation rate was applied for population-level estimates for two reasons: first, the divergence rate with 13.4%/Myr is near the high end of estimates reported for sardines fishes (Bowen et al., 2006); second, Ho et al. (2008) reported up to 10 × elevated rates on short time scales, so an order-of-magnitude range of rate estimates should bracket the most extreme plausible cases. We used the AMOVA analysis to estimate the population structure of C. nasus throughout the distribution of species (Excoffier et al., 1992). According to the geographic pattern, three groups, including the middle reaches of the Yangtze River group, the low reaches of the Yangtze River and the Chinese coastal group, and the Ariake Sea group, were designed for the AMOVA analyses. Ten thousands permutations were used to test the significance of covariances component from different assembly levels. Pairwise different values among the populations were also carried out in the examined ranges (Excoffier et al., 1992). Ten thousands permutations were also used to test the significance of the ST for each pairwise comparison. All the parameters of the population structure were calculated using the ARLEQUIN software.
We used the Google Earth (Version 6) to obtain the accurate geographical distances along with coast lines among localities. Then, we used Mantel tests implemented in GeneAlex to calculate the corelations between ST /(1-ST ) and standardized geographical distance (Peakall and Smouse, 2006). We used the neutrality test, mismatch distribution analyses implemented in the ARLEQUIN to check population historical demography for C. nasus based on the mtDNA D-loop (Fu, 1997). The Bayesian skyline analyses were used to estimate the fluctuations of effective population sizes (ESS) for C. nasus . We used the beast 1.7.5 software to create the Bayesian skyline plot (BSP) based on the mtDNA D-loop datasets (Drummond et al., 2012). Genealogies were combined from 1 × 10 8 steps with a burn-in of 1 × 10 7 lineages. To obtain the effective convergence, we performed multiple analyses that were run for 5 × 10 8 iterations with a burn-in of 107 under the HKY + I + G model, a strict molecular clock and stepwise skyline model. The results of 10 times runs for C. nasus were integrated using LogCombiner software, and the values of ESS for each parameter all exceeded 200. The skyline plot was generated by Tracer 1.5 (Rambaut and Drummond, 2007).

Genetic Diversity of C. nasus
The mtDNA D-loop (624-703 bp) of C. nasus contained one to four copies of tandemly repeated 38 bp sequence units. After the removal of the tandemly repeated sequences and a conserved 10 bp portion of the tRNA Pro gene, nucleotide sequence with length of 465 bp mtDNA D-loop was determined in 424 individuals. The mtDNA D-loop alignment consisted of 108 polymorphic sites, of which 88 were transition sites, 19 were transversion sites, and 7 were indels, yielding 171 haplotypes. One hundred and twenty-four haplotypes (72.5%) were represented by a single sequence. Forty-one of the remaining haplotypes were shared by more than one population and the last remaining 6 haplotypes were shared by at least two individuals in the same population. The average values of intrapopulation genetic diversity indices showed that a pattern with both high levels of haplotype diversity (h = 0.97) and nucleotide diversity (π = 0.02) were detected in the examined range.
The mtDNA Cyt b data set, with an aligned fragment length of 402 bp, was obtained from 139 individual from 8 localities. A total of 9 substitutions were detected in the mtDNA Cyt b fragment with 8 transitions and 1 transversions, which defined 10 haplotypes. Seven of the ten haplotypes were represented by a single individual. The average values of haplotype diversity (h) and nucleotide diversity (π) for all the examined populations were 0.455 and 0.001, respectively.

Phylogenetic Relationships and Four Lineages
Four genealogical clades (labeled A, B, C, and D) were detected in the phylogenetic tree under the HKY + I + G model based on the mtDNA D-loop data set (Figure 2). The MST analysis also indicated that four distinct lineages with three or six mutational steps corresponding to the phylogenetic tree were existed in the examined range (Figure 3). The MST was characterized by star-like structure with common haplotypes distributing in the center of the structure (Figure 3). The genetic distance among four clades were also calculated as follows: A/B, 0.0270; A/C, 0.0236; A/D, 0.0296; B/C, 0.0145; B/D, 0.0252; and C/D, 0.0220. According to the population-level divergence rate (13.4% per Myr), the divergence events among these genealogical clades happened at about 110,000-220,000 years ago.
The distribution of haplotype frequencies for the four genealogical clades showed significantly geographic differentiation in the examined ranges. The populations along the middle reaches of the Yangtze River were mainly occupied by the lineage A, which included 44 haplotypes with 132 individuals. In addition, low frequencies of the haplotypes from lineage A were also detected in the six populations from the lower  reaches of the Yangtze River and the Huai River (Table 1 and Figure 1). Additionally, two haplotypes of lineages A (hap97 and hap98) were also found in Dandong and Wenzhou. Five mutational steps existed between these two haplotypes and other haplotypes of lineage A (Figure 3). The frequency of lineage B (29 haplotypes, 41 individuals), which was sympatric with lineage D (92 haplotypes, 237 individuals) in the populations from the coast of China, declined along the coast from Dandong to Wenzhou. The lineage B haplotypes were also observed at low frequency in those populations from the middle and lower reaches of the Yangtze River and the Huai River system, such as Poyang Lake, Chao Lake, Nanjing, Nantong, and Zhoukou. Only lineage B haplotypes were detected in the Ariake Bay. The frequency of lineage D was higher than that of lineage B in the examined populations. Lineage D also was found at low frequency in populations from the middle reaches of the Yangtze River. The haplotypes of lineage C (6 haplotypes, 14 individuals) were distributed in populations from the coast of China except Dandong, and with a high frequency (33.3%) in Dongying. The lineage C haplotypes also were found at a low frequency in the Yangtze River basin, such as in Dongting Lake, Nanjing and Zhenjiang.
For the mtDNA Cyt b sequences, the MST of haplotypes revealed the C. nasus and C. grayii were separated from each other by 10 mutational steps (d = 2.5%). The MST for mtDNA Cyt b haplotypes was also characterized by star-like structure for C. nasus and C. grayii (Figure 4). For C. nasus, two main haplotypes (hap1 and hap3) were separated from each other only by one mutational step (Figure 4). For C. grayii, 8 haplotypes were also divergent from one to another by one mutational step. Furthermore, the MST of C. grayii also suggested the sequences of C. nasus from Ma et al. (2010) and Qiao et al. (2013) were within the lineage of C. grayii (Figure 4). As a consequence, these sequences reported by Ma et al. (2010) and Qiao et al. (2013) should belong to the specimen of C. grayii.

Population Structure
AMOVAs based on the geographic scale were carried out and revealed 50.92% variation among groups (P < 0.01) and 1.32% differentiation among populations (P < 0.01). The values of CT and SC were calculated as follows: 0.558 (P < 0.01) and 0.022 (P < 0.01). The pairwise ST values showed that non-significant differentiations were checked among populations from the middle reaches of the Yangtze River group (Table 3). Within the lower reaches of the Yangtze River and the Chinese coastal group, the genetic differentiation was significant between Dandong and other most of the populations from the Huai River, the lower reaches of the Yangtze and Wenzhou. The genetic A significant low-level signal of isolation-by-distance (IBD) was observed among all populations (R 2 = 0.2175; P = 0.001; Figure 5A). The significant low-level IBD pattern (R 2 = 0.1978; P = 0.018) was also found among all populations except the five populations of the middle reaches of the Yangtze River group (R 2 = 0.0161; P = 0.338; Figures 5B,C). No IBD pattern was detected for the twelve populations from the lower reaches of the Yangtze and the coast of China group (R 2 = 0.056, P = 0.123; Figure 5D).

Population Demography
The F S tests for all the samples and the four lineages were significantly negative (Table 4), indicating a possible population expansion event for the C. nasus. The unimodal mismatch distributions were detected in each lineage A, B, and C, which further demonstrated that the population expansion event had happened in the species (Figure 6). The bimodal mismatch distribution was detected in the lineage D, one of which corresponded to the differences among subclades, and the other of which corresponded to the differences among individuals within subclades (Figure 6). Based on τ values and 10-fold faster rate (13.4% per Myr), the expansion time for the three clades was dated back to the late Pleistocene period. Bayesian skyline plots demonstrated demographic expansion in each of the four lineages (Figure 7). Applying the 10 × faster mutation rate (13.4%/Myr), the population expansion time of lineage A occurred at about 20,000-30,000 years ago. For lineages C and D, the demographic expansion was approximately 20,000-25,000 and 25,000-27,000 years ago, respectively. Furthermore, the expansion time for lineage B was much higher than those of the left lineages happened at about 55,000-75,000 years ago.

DISCUSSION
The fossil records were always used as calibration points in the determination of the divergence rate, paleogeologic range, and paleogeological environment of the organisms (Lavoué et al., 2012(Lavoué et al., , 2013. Lavoué et al. (2013) investigated 73 clupeoid taxa collected from 5 families and used fossils as the calibration to estimate the node age. According to their analyses, we could identify the divergence time between the C. nasus and C. ectenes was as approximately 3.6 Myr. The lineage-specific divergence rate for mtDNA D-loop sequence was approximately 1.34%/Myr for the Japanese grenadier anchovies based on the comparison results of mtDNA Cty b.
Recently, Ho et al. (2005) argued that the phylogenetic substitution rates could not be directly applied to populationlevel analyses. The substitution rates should be ensured based on the calibration point age (Emerson, 2007;Burridge et al., 2008). To account to some extent for such concerns, Liu et al. (2011b) applied two different alternative rate calibrations ("phylogenetic" and 10-fold faster rates) to the analyses of Pleistocene lineage diversification for Pacific herring. However, the substitution rates for mtDNA D-loop was uncertain and seemed vary among taxonomic groups (Bowen et al., 2006). Various molecular rates were used in the previous researches of Clupeoidei as follows: 7.2%/MY for California anchovy (Engraulidae), 5-10%/MY for Japanese anchovy and Australian anchovy (Engraulidae), and 15-20%/MY for Indo-Pacific sardines (Clupeidae; Bowen and Grant, 1997;Liu et al., 2006;Díaz-Viloria et al., 2012). Thus, the "phylogenetic" rate (1.34%/Myr) based on the presumed split between C. nasus and C. grayii was slower and inappropriate for population-level analyses of C. nasus. In this study, we applied a substitution rate (13.4%/Myr) that is 10-fold higher than the phylogenetic rate and is close to the estimates reported for sardines (Bowen et al., 2006). Four major genealogical clades were detected in C. nasus based on the mtDNA D-loop analyses. The freshwater form of C. nasus mainly existed in the Yangtze River under the effects of geographic isolation, which strongly supported by the geographical distribution of lineage A. The lineage B of C. nasus was originated from the geographic isolation between Ariake Bay and Northwestern Pacific during the late Pleistocene. The lineage D of C. nasus was mainly composed of populations from the Chinese coast waters and the lower reaches of the Yangtze River. Lineage C may be constituted by some ancient haplotypes remained in the East China Sea following the last glacial in accordance with a low frequency distribution along the Chinese coast populations and the low reaches of the Yangtze River. According to the population-level divergence rates utilized in the present study, lineages A, B, C and D separated from one another at approximately 110,000-220,000 years ago.
Additionally, two dominant haplotypes were found in analysis of the partial mtDNA cytochrome b sequences. One was mainly shared by the five populations in the middle reaches of the Yangtze River and the other was shared by the populations along coast and the lower reaches of the Yangtze River. The results of phylogenetic analysis for both the mtDNA D-loop and Cyt b sequences were consistent with the allocation of the freshwater form and anadromous form reported by Cheng (2011).
The genetic variation of marine species has been shaped by the periodic climatic oscillations of the Pleistocene in East Asia, which have raised range contractions, secondary contact between different basins . Strong genetic divergence between populations in the middle and lower reaches of the Yangtze River have also been reported for Macrobrachium nipponense, Neosalanx taihuensis, Hypophthalmichthys molitrix, Ctenopharyngodon piceus, Aristichthys nobilis, and Mylopharyngodon piceus (Lu et al., 1997;Li and Lu, 1998;Feng et al., 2008;Zhao et al., 2008). These concordant results showed that paleogeographic factors had greatly influenced the evolutionary processes of species in the middle and lower reaches of the Yangtze River. The Pleistocene glacial-interglacial cycles was characterized by the intervals of ∼100 kyr over the past ∼800 kyr (Lambeck et al., 2002). During the glaciations periods, the sea level reduced to the 120-140 m and the basin in the middle of the Yangtze River formed deeply incised river valleys (Faure et al., 2002). During the inter-glaciations periods, water level of the Yangtze River rose and the lakes expanded after heavy rainfall (Yang, 1986). The lower reaches of the Yangtze River could be inundated by the seawater during the sea level rise.
The Yangtze River was an inland river derived from the impacts of arid climate, reduced rainfall and desertification during the LGM, which isolated with East China Sea (Xiao et al., 2003). Thus, some populations of C. nasus inhabited along in the deep-water channel of the Yangtze River during the glaciation periods and became the freshwater form under the effect of geographic isolation. The glacial refugium for the freshwater  population was inferred to be located in the middle reaches of the Yangtze River. However, the existence of special sediment proves that Dongting Lake and Poyang Lake were formed less than 10,000 years ago (Guan and Cai, 1986). Therefore, Dongting Lake and Poyang Lake were unlikely to have been the glacial refuges for these freshwater populations. The deficiency of gene flow between the migratory and non-migratory populations caused the differentiations of morphology, behavior, and genetics, and further gave rise to incipient speciation of island freshwater fish (Michel et al., 2008). Freshwater resident individuals of C. nasus have short maxilla and small body size, and their main breeding season is different from the individuals in the lower reaches of the Yangtze River (Yuan, 1986). The anadromous populations in the lower reaches of the Yangtze River mainly spawn from April to May, whereas the freshwater populations in the middle reaches of the Yangtze River spawns from May to August (Zhang et al., 2005;Luo, 2006;Liu, 2008). The incongruence of the breeding season probably accelerates the genetic divergence between the freshwater form and the anadromous form. Strong genetic differentiations between the Ariake Sea populations and the Chinese coast populations have been observed in some marine fishes, such as Pennahia argentata, Verasper variegatus, Salanx ariakensis, and Mugil cephalus, marine invertebrates and marine mammal (Han et al., 2008;Yang et al., 2008;Hua et al., 2009;Mao et al., 2011;Sekino et al., 2011;Shen et al., 2011). The Ariake Bay is a semi-closed basin, whose geographical environment has been greatly influenced by the sea-level fluctuations during the Pleistocene glacial-interglacial periods. The population genetic compositions of marine species in the Ariake Bay have been proved to be isolated with the other populations distributed along the extension of the bay (Shimoyama et al., 1996;Kikuchi, 2000;Sato, 2001;Kojima et al., 2004). In the present study, the results of the distribution of lineages B and D strongly supported the hypothesis that the population genetic composition of C. nasus between Chinese coast and Ariake Bay of Japan had been dramatically shaped by the Pleistocene climatic oscillations. The results of lineages B and C showed deeply separated mtDNA lineages coexisted in the examined ranges, however, no significantly geographical isolation could be responsible for the sympatric differentiations.
Although it was reported that enough population size with stable or grew throughout long history would accumulate numerous sequence differences (Avise, 2004), we suggested that the secondary contacts of different lineages originated from the isolations among different glacial-refugia (e.g., the Ariake Bay, the East China Sea, the Huang River, and the Huai River) during the Pleistocene period should be responsible for the geographic distribution of the lineages in C. nasus.
The current analyses do not support the interspecific genetic divergence as reported by Ma et al. (2010). Their conclusions suggested that there were some highly diverged mtDNA cytochrome b haplotypes (29 substitutions) in Fuzhou (FZ) and Wuhan (WH) populations. But from our analyses, we can conclude the special haplotypes in Fuzhou and Wuhan populations belong to the specimen of C. grayii, not C. nasus (Ma et al., 2010). Gray's grenadier anchovy, C. grayii are widely distributed in the East and South China Sea (from Fuzhou City to Guangxi Province; Whitehead et al., 1988;Zhang, 2001). The Gray's grenadier anchovy doesn't inhabit along the reaches of Yangtze River. Thus, Ma et al. misidentified the specimen and confused the sample sites of C. nasus and C. grayii.
Strong gene frequency changes of the three dominant lineages (A, B, and D) in the examined ranges showed highly limited genetic exchange existed among groups of C. nasus. Our results also showed that the presence of lineage A (freshwater individuals) occupied low frequencies of haplotypes in the low reaches of the Yangtze implied extremely limited gene flow existed between the anadromous group and freshwater group. Our result was also supported by the phylogenetic analysis of representative individuals clearly corresponding to the anadromous group and freshwater group of C. nasus based on the entire genome SNP set (Xu et al., 2020). The middle reaches of the Yangtze River harbors a large number of lakes with different sizes (e.g., Tianezhou, Poyang Lake, and Dongting Lake), which would support optimum living and spawning conditions (e.g., fresh water, water temperature (15∼27.5 • C), depth (0.5∼3.0 m), suitable current (0.057∼0.075 m/s), and abundant zooplankton required for spawning) for freshwater group of C. nasus. The optimum living and spawning conditions for freshwater group of C. nasus, especially high quality fresh water, gentle currents and abundant zooplankton, could not be found in the low reaches of the Yangtze River and the China coast (Xu et al., 2020). The effect of isolation (the middle reaches of the Yangtze River vs. the low reaches of the Yangtze River and the China coast) and adaptation to the new freshwater habitats of the middle reaches of the Yangtze River for C. nasus freshwater form might be responsible for the module of finitely asymmetric gene flow from freshwater to anadromous populations (Zhang, 2001;Schluter and Conte, 2009). The conclusions were also further supported by the module with high gene exchanges among anadromous form populations for C. nasus (Tonteri et al., 2007;Bourret et al., 2013a,b;Perrier et al., 2013). The lineage A was also found in the populations (ZK and HL) from the system of the Huai River, which might be derived from the gene flow channel (Hongze Lake) connecting the Huai River and the Yangtze River since 1851 (Lu and Wang, 2009). Therefore, individuals of lineage A might have dispersed into the system of the Huai River by the connection between the two rivers. Then, the freshwater form (lineage A) was amplified in the Huai River to bring about the high frequency of lineage A existed in the Zhoukou population, which can be explained by the "transporter" hypothesis. Moreover, the presence of lineages B, C, and D at low frequencies within the five freshwater populations may suggest hybridization between the two forms, or upstream spawning migration of anadromous individuals. High frequencies variations of haplotypes were detected in the lineages B and D between the Ariake Bay and the coast of China, which suggested highly limited gene flow existed among populations. During the Pleistocene ice ages, the East China Sea was pushed to the East to be a huge enclosed sea because of exposed continental shelves (Liu et al., 2007). Therefore, the populations of C. nasus in the Ariake Sea were separated from those populations in the East China Sea. The isolation of marginal seas and the contact and mix among populations during the Pleistocene glacial and interglacial epoch would be responsible for the current distribution of lineages B, C and D in the East China Sea.
All the ST values of the comparisons among five populations of the middle reaches of the Yangtze River group were not statistically significant, which indicated the genetic structures among the freshwater forms were lacked. The result was further supported by the IBD pattern (R 2 = 0.0161, P > 0.05) in the middle reaches of the Yangtze River ( Figure 5B). A limited range of activities and current coupled with pelagic eggs of C. nasus could be responsible for the frequent gene flow among populations in the middle reaches of the Yangtze River (Yuan et al., 1976;Ni et al., 1990;Xu et al., 2020). The similar results with no significant IBD pattern was also detected among populations from the lower reaches of the Yangtze and the coast of China group (R 2 = 0.056, P = 0.123), which might be derived from integrations of the upstream spawning migration behavior of anadromous individuals and transportation of juveniles by China Coastal Current (Figure 5D; Yuan et al., 1976;Yuan et al., 1980;Liu et al., 2007;Cheng, 2011;Xu et al., 2020). The Hongze Lake and Zhoukou populations belonged to the Huai River system which was connected by the Yangtze River via a human-made channel. However, the dams located the humanmade channel could limit the Yangtze River populations to enter into the Huai River. Thus, the existence of genetic substructure between the Huai River and the Yangtze system was supported by the significant ST values of three comparisons between populations (HL & NJ, ZK & NJ, and ZK & TL). Significant ST values were observed between the Dandong (the estuary of the Yalu River) population and other three populations (NT, TL, and CL) from the lower reaches of the Yangtze River. This level of subdivision may be reflecting the homing behavior of adults (Campton and Utter, 1987;McDowall, 2001;Keefer and Caudill, 2014). Furthermore, the significant ST  Therefore, although no IBD pattern was detected for populations from the lower reaches of the Yangtze and the coast of China group, different homing behavioral mechanisms may bring about the weak genetic structure between the ZK population and other two estuary populations (DY and PJ).
Neutrality test, mismatch distribution and network analyses combined with Bayesian skyline plots consistently showed that all the four lineages (lineages A, B, C, and D) had populations with relatively stable over time and then subjected to sudden expansion during the late-Pleistocene era (c. 20,000-75,000 years ago). According to the divergence rate of 13.4%/Myr, the initial population expansion of lineages A, C and D occurred 20,000-30,000 years ago that largely coincided with the onset of the current Pleistocene interglacial period. During the last interglacial period, the transgression has happened in the East China Sea Shelf, which caused the coastline to migrate from the Okinawa Trough to the coast of the Bohai bay by 1200 km (Wang, 1999;Xu and Oda, 1999). Under this scenario, populations of C. nasus from the China coast would expand their range along the coastal areas following the glacial retreat. Moreover, water level of the Yangtze River rose and the lakes expanded after heavy rainfall during the last deglaciation. Freshwater populations in the Yangtze River might also expand their range.
Our results showed that the phylogeographic and demographic history of C. nasus has been greatly influenced by the Pleistocene glacial-interglacial cycles. With respect to phylogeography, the genetic divergences between four lineages were observed in the middle reaches of Yangtze River, Chinese coastal and Japanese coastal populations. Geographic isolation resulted from climatic fluctuations of the Pleistocene has a great influence on the distribution of the four lineages. The glacial refugium for the freshwater population of C. nasus was proved to be located in the middle reaches of the Yangtze River. With respect to the historical demography, our result showed that C. nasus had experienced dramatic population expansion event during the glacial period of the Pleistocene. But beyond all that, additional nuclear markers are needed to be used for further studies to test our conclusion of the biogeographical pattern for C. nasus.

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 Institutional Ethics Committee of Center for Ocean Mega-Science.

AUTHOR CONTRIBUTIONS
TG, QY, and YX conceived and designed research. TG, QY, NS, YY, and YX conducted experiments, analyzed data, and wrote the manuscript. Authors critically reviewed and approved the manuscript. All authors contributed to the article and approved the submitted version.