1 INTRODUCTION

The Japanese eel (Anguilla japonica) is a temperate catadromous anguillid eel, distributed in the rivers of China (China's mainland and Taiwan region), Japan, and Korea (Tesch, 1977). It spawns in the waters to the west of Mariana Islands near 14°N–16°N, 142°E (Tsukamoto, 1992, 2006), 2 000 to 3 500 km away from the East Asian continent. The leptocephali hatch between April and November (Tsukamoto, 1990; Tzeng, 1990; Tsukamoto et al., 2003), and then drift with the currents, reaching the coasts of East Asia in 4 to 6 months. They metamorphose into glass eels along the continental shelf and then enter estuaries where they continue to grow as elvers (Tesch, 2003). Elvers grow to yellow eels during their upstream migration, and yellow eels spend the next 5 to 10 years in freshwater rivers. Once eels attain sexual maturity, they migrate back to spawning ground for reproduction (Tesch, 1977).

The Japanese eel is a commercially important fish species in East Asia. And the stocks are outside safe biological limits and the fishery has not been sustainable in recent years. The Japanese eel population is currently estimated to be less than 10% of 1970s level (Dekker, 2003), and it is classified as "Endangered" by the International Union for Conservation of Nature (IUCN) Red List (https://www.iucnredlist.org/species/166184/1117791). Although the exact reason for the decline is unclear, there is a direct correlation between the decline in eel catch and population size to the habitat modifications and impediments to upstream migration (Itakura et al., 2015). Other reasons such as pollution, disease, oceanic climate change and habitat losses could also influence the abundance of the Japanese eel population (Arai, 2014). However, the reliance of eel aquaculture on wild-caught juveniles or glass eel seems to be more important reason (Tsukamoto, 2013). Due to the large distribution range of A.japonica, the effective management needs to involve the efforts from all countries within its range. And clarifying the genetic population structure of A.japonica is the basis of multinational cooperation on management and protection. In addition, it is also of vital importance to understand the genetic consequences of external disturbance to this endangered species, such as on-going climate change, ocean warming and overfishing.

Some work has been conducted to test panmixia of A. japonica using different molecular markers; however, the results are not consistent with each other. Panmixia was accepted based on analyses on mtDNA sequences (Sang et al., 1994; Ishikawa et al., 2001). Microsatellite markers, which have more power in studying population structure, were also used to solve this problem. Tseng et al. (2006) divided the populations of A. japonica into two groups based on eight microsatellite markers, and they found no isolation by distance (IBD) pattern in either group. The results challenged the panmixia hypothesis by finding small but significant genetic differences. However, the study used a small number of markers and reported higher temporal than geographic difference. Notwithstanding the results have rejected the null hypothesis of panmixia for the genetic structure of anguillid eels, several studies have challenged this opinion recently. Han et al. (2010) used eight microsatellite markers analyzing the spatial and temporal structure of A. japonica. Although a significant difference was seen among annual cohorts and spatial samples, no specific temporal or spatial scale patterns were observed in the pairwise F_ST tests or the phylogenetic tree of all samples and the results of the IBD test and the isolation by time (IBT) test were both insignificant. A stable genetic structure could not be observed. The result provided evidence for panmixia. After that, Minegishi et al. (2012) also found no significant genetic difference among 9 localities using 6 microsatellite markers and reported that Japanese forms a panmictic population.

Traditional approach using neutral markers have reported shallow population structure in many marine fishes, which is assumed to be the result of high levels of gene flow. However, in species with large effective population size, the weak genetic structure could also result from the limited effects of genetic drift (Hauser and Carvalho, 2008). Russello et al. (2012) detected ecotype-level divergence using 8 outlier loci in Okanagan Lake kokanee (Oncorhynchus nerka); however, there was no evidence of divergence at neutral loci. Milano et al. (2014) detected a dramatic divergence between Atlantic and Mediterranean populations and fine-scale significant population structure in European hake (Merluccius merluccius) using outlier loci, which was not revealed using neutral loci. Many marine organisms have relatively large effective population size, and they are under limited effects of genetic drift. In such a scenario, gene-associated markers might have more power in unveiling the cryptic genetic population structure than neutral markers.

Local adaptation is how organism response to selective pressures in their local habitats, acting on genetically controlled fitness differences among individuals (Rellstab et al., 2015). Since eels are panmictic and show no philopatry (Als et al., 2011; Pujolar et al., 2014), long-term local adaptation is not possible in eels but single-generation signatures of local selection still can be detected (Gagnaire et al., 2012; Pujolar et al., 2014), which is a completely different scenario relative to other species. In other fish species, fish can adapt to local conditions and those individuals more resistant to adverse environmental condition in a population will survive, while less adapted individuals will not, and this is passed on to the next generation. However, this is not the case in eels: Individuals may have some good genetic combinations which provide good adaptation for survival in their freshwater habitat, however, habitats for offspring of these individuals would change due to random dispersal, which could lead to the loss of these adaptations. Trans-generational signatures of selection cannot be detected while within-generation selection acting on genes could be found. Recently, gene-associated markers have shown their power in detecting signature of selection in several organisms. In the case of anguillid eels, signature of selection is much more difficult to find, due to the special migratory features. However, recent studies have found signature of selection in North Atlantic eels. Gagnaire et al. (2012) scanned for footprints of selection in the American eel (Anguilla rostrata) using a panel of 100 single nucleotide polymorphisms (SNPs) and found the allele frequencies of thirteen loci had a correlation with environmental factors. In the study of European eel (Anguilla anguilla) using the same SNP panel, 10 loci were detected showing genetic-by-environment associations (Ulrik et al., 2014). Pujolar et al. (2014) detected signature of selection in the European eel using 50 354 restriction-site-associated DNA sequencing (RAD) generated SNPs, and dozens of loci showed evidence of divergent selection. However, the study on local adaptation is still absent in A. japonica.

The previous studies have led to conflicting results on panmixia, but all previous studies used a very limited (6–8) number of microsatellites. In the present study, 24 polymorphic gene-associated microsatellite markers were developed for A. japonica. These markers were then used to test the hypothesis of panmixia and to detect footprints of selection in A. japonica.

2 MATERIAL AND METHOD 2.1 Sampling and DNA extraction

Glass eels or yellow eels were collected from 13 sampling locations spanning the natural species distribution between January 2011 and April 2016 (Table 1), and the whole organism of these eels was preserved in 95% ethanol. Whole genomic DNA was extracted using a standard phenol-chloroform extraction and ethanol precipitation, following the proteinase K digestion.

2.2 Microsatellite markers development and genotyping

Microsatellite loci were obtained in two ways. Firstly, we searched for the transcript sequences of the genes against the online transcriptome database for A. japonica (http://molas.iis.sinica.edu.tw/jpeel/) built by Hsu et al. (2015) using the names of selected genes in A. rostrata (Gagnaire et al., 2012) and A. anguilla (Pujolar et al., 2014) as keyword. The detail information of the microsatellite loci in each transcript sequence were obtained using MISA software (http://pgrc.ipk-gatersleben.de/misa). Secondly, we developed another panel of microsatellite loci from the transcript sequences randomly chosen from the database of A. japonica. To obtain a longer flanking region for primer design, each transcript sequence with the microsatellite locus was mapped on the draft genome of A. japonica (Henkel et al., 2012). Finally, Premier 5.0 software (Premier Biosoft International, Palo Alto, CA, USA) was used to design primers. 23 primer pairs were designed on the basis of the gene under selection in North Atlantic eels and 12 primer pairs were designed from transcript sequence randomly selected from database. Out of the 36 primer sets, 24 primer pairs consistently amplified fragments of the expected size and were polymorphic. Among the 24 polymorphic loci, 18 loci were developed associated with the gene under selection in North Atlantic eels, and the other 6 loci were developed based on transcript sequence which were randomly selected from the database. PCR was carried out using a three-primer system (Liu and Avise, 2011). All forward primers were 5′ tailed with a M13 forward sequence (5′-GGAAACAGCTATGACCAT-3′). An M13 primer of the same sequence but 5′-labelled with three different fluorescent dyes was used in combination with forward primers. All loci were amplified separately on a Mastercycler (Eppendorf). The PCR consisted of approximately 50 ng genomic DNA, 2 μL of 5× buffer (Promega, Madison, WI), 0.2 mmol/L of each dNTP, 0.2 μmol/L labeled M13 reverse primer and locus specific primer without tail, 0.02 μmol/L locus specific primer with M13 reverse tail, 0.25 U Taq DNA polymerase (Promega, Madison, WI) and water up to 10 μL. The thermal cycling protocol consisted of 3 min at 95℃ followed by 35 cycles of 20 s at 95℃, 52 s at 52℃ and 30 s at 72℃, and a final elongation of 10 min at 72℃. PCR products labelled with different fluorescent primers were pooled and electrophoresis was made on an ABI 3730xl automated sequencer (Applied Biosystems). Size determination of alleles were made by comparison with the GS500-ROX size standard using GeneMarker (SoftGenetics, State College, USA).

2.3 Genetic diversity analysis

Genetic diversity indices including observed heterozygosity (H_O), expected heterozygosity (H_E), the number of alleles (N_a), and inbreeding coefficient (F_IS) were calculated using GenAlEX 6.5 (Peakall and Smouse, 2006, 2012). The presence of null alleles was assessed using MICRO-CHECKER software version 2.2.3 (Van Oosterhout et al., 2004). Genotypic linkage disequilibrium (LD) and deviation from Hardy-Weinberg equilibrium (HWE) were estimated using GENEPOP 4.5.1 (Rousset, 2008). The significant thresholds for HWE and LD tests were adjusted using Bonferroni correction.

2.4 Population structure

Global F_ST value across all populations was calculated using FSTAT version 2.9.3.2 (Goudet, 2001). Pairwise F_ST values between all populations were obtained using ARLEQUIN 3.5 (Excoffier and Lischer, 2010) and significance was assessed by 10 000 bootstrapping permutations. A modified false discovery rate (FDR) correction was applied to adjust the significant threshold of pairwise F_ST values (Benjamini and Yekutieli, 2001). The Bayesian clustering program STRUCTURE version 2.3.4 (Pritchard et al., 2000), was used to infer population structure and assign individuals to populations. Admixture model was used, assuming correlated allele frequency, with the option of 'with no prior knowledge of sampling locations'. The program was run with an initial burn-in of 100 000 cycles and 1 000 000 Markov chain Monte Carlo (MCMC) repetitions, and K was modeled from K=1 to K=13. Runs were iterated 10 times for each K. The optimal value of ΔK was determined according to Evanno et al. (2005) method using StructureSelector (Li and Liu, 2018).

2.5 Footprints of selection

The F_ST outlier approach (Beaumont and Nichols, 1996) implemented in LOSITAN (Antao et al., 2008) was used to identify outlier loci. The program generated a plot of F_ST vs heterozygosity.

One hundred thousand simulations were performed using the stepwise mutation model (SMM) and the neutrality of microsatellite markers was determined. Markers having F_ST values higher than the 95% confidence interval were inferred to be subject to divergent selection, and markers having F_ST values lower than the 95% confidence interval were inferred to be subject to balancing selection. False discovery rate was set to 0.01.

Spatial analysis method (SAM) was used to assess the association between allelic frequencies and environmental variables, based on multiple univariate logistic regression (Joost et al., 2007). Environmental factors were chosen including degrees north latitude, degrees east longitude, and sea-surface temperature at river mouth averaged across the 10 days, 30 days and 3 months that preceded the sampling date. Temperature data were acquired from a National Oceanic and Atmospheric Administration (NOAA) database containing geo-referenced sea-surface temperatures. The FY population was excluded from this analysis, because the sample included different development stages and the sampling date was also different. Two populations (MJ and MC) from Japan and one population from Taiwan, China (MY) were excluded from genotype-temperature association analysis, because the samples were collected from different date in the same month.

3 RESULT 3.1 Genetic diversity

The average number of alleles per locus across populations ranged from 4.846 (SE=0.274) to 29.231 (SE=1.156) (Table 2). The H_O and H_E of each sample ranged from 0.675 (SE=0.038) and 0.746 (SE=0.039) to 0.716 (SE=0.033) and 0.769 (SE=0.034) (Table 3). Among a total of 312 sample x locus combinations, 29 were deviated from HWE after Bonferroni corrections (α=0.05, K=13). After Bonferroni correction, no linkage disequilibrium was detected (α=0.05, K=276). Null alleles existed at eight loci in more than three samples (not in all samples).

3.2 Population structure

The global non-significant F_ST value of -0.001 (95% CI: -0.002 < F_ST < 0.000) indicated that there was no significant genetic differentiation among these samples. None of the pairwise F_ST values were significant after FDR correction (α=0.05, K=78) (Table 4). STRUCTURE analysis identified K=6 was the optimal number of clusters (Fig. 1). However, all populations showed similar patterns when K=6, demonstrating no genetic structure (Fig. 2).

3.3 Test of selection

No outlier was found in outlier analysis (Fig. 3). No significant correlations were found in SAM, indicating that the 24 loci were not under spatially varying selection.

4 DISCUSSION 4.1 Confirmed panmixia in A. japonica

Japanese eel has long been considered having a panmictic population structure because sexually matured eels from their distribution range migrate and spawn in a single spawning site and the larvae disperse from the spawning site via the North Equatorial Current (NEC) and the Kuroshio Current (KC) and are transported back to the freshwater habitats randomly. However, there is still debate on the hypothesis of panmixia for A. japonica. Chan et al. (1997) detected spatial genetic heterogeneity in A. japonica using allozyme markers, the significant clinal shift in allelic frequency was putatively attributed to single-generation selection along an environment gradient (Avise, 2003). The geographic range of the sampling sites in our study was similar to that in Chan et al. (1997), however, we found no spatial heterogeneity in the gene frequency of the 24 microsatellites. Tseng et al. (2006) divided the populations of A. japonica into two groups including a low-latitude group (South China and Taiwan, China) and a high-latitude group (Japan, Korea, and North China) using eight microsatellite markers. However, their study used a small number of markers and reported higher temporal difference than geographic difference. Han et al. (2010) suggested that the genetic partitioning detected in Tseng et al. (2006) could be the results of 'chaotic genetic patchiness' and random variations in parental contributions to reproductive activity, incomplete mixing of larvae and kin aggregation might contribute to the subdivision of the Japanese eel population. Han et al. (2010) conducted the most comprehensive study on A. japonica using microsatellite markers and confirmed panmixia of this species.

In the present study, all measures of differentiation were accordant with a panmixia scenario. Our results supported the hypothesis of panmixia in A. japonica. Great efforts have been put into the study of genetic population structure for A. japonica, however, no one has ever tried to use gene-associated markers. Gene-associated makers can be a good option to detect subtle genetic population differentiation in marine organisms. On one hand, traditional neutral markers may be ineffective when populations are recently isolated, and divergence is not yet reflected at neutral loci (Russello et al., 2012). On the other hand, genetic drift has little effects in species with large effective population size (such as marine fishes), which may weaken the power of neutral markers (Milano et al., 2014). At gene-associated markers putatively under selection, genetic differentiation could be larger than that at loci in which among-population differences are caused by genetic drift (Holsinger and Weir, 2009).

There are several studies showing the power of outlier loci in detecting population structure. Russello et al. (2012) successfully detected ecotype-level divergence using eight outlier microsatellite markers, whereas there was no evidence of divergence at neutral loci. By using outlier loci, Milano et al. (2014) revealed a dramatic divergence between Atlantic and Mediterranean populations of European hake and fine-scale significant population structure among regional populations, which was not detected by using neutral markers. Another advantage of the markers in this study is their potential transferability across taxonomic boundaries. Compared with traditional microsatellite markers, microsatellite markers based on transcript sequences are potentially more transferable across taxonomic boundaries (Chagné et al., 2004; Liewlaksaneeyanawin et al., 2004; Gutierrez et al., 2005; Pashley et al., 2006).

For the first time, our study confirmed the panmixia of A. japonica using gene-associated microsatellite markers. We suggested that A. japonica should be managed as a single unit. Management and conservation efforts must be coordinated at the transnational level, as over-exploitation in any region will decrease recruitment across its whole distributional range.

4.2 Different pattern of spatially varying selection from North Atlantic eels

Both LOSITAN and SAM analysis detected no signals of spatially varying selection. These results indicated that the 24 loci were not targets of spatially varying selection and were nearly neutral genetic markers for Japanese eel. Since most of the microsatellite markers used are associated with the genes under selection in the North Atlantic eels, our results suggested that different selection pattern existed between Japanese eel and the North Atlantic eels. Furthermore, two recent studies found no overlap of loci under selection in the two North Atlantic eel species either. Gagnaire et al. (2012) scanned for footprints of selection in the American eel using a panel of 100 SNPs and 13 loci showed significant correlations between allele frequencies and environmental variables. Ulrik et al. (2014) found evidence for spatially varying selection at 10 loci in the European eel using the same SNP panel; however, none of these were the same loci that showed significant associations in the American eel. The contrasting results in the two species suggested that there were two different patterns of selection in the North Atlantic eels, at least at the level of the individual genes assessed (Hemmer-Hansen et al., 2014). Our results suggested that A. japonica might occupy a third selection pattern, which was different from both the North Atlantic eel species. Due to the latitudinal (rather than longitudinal) spread of freshwater habitats and the longer migration distance to the spawning site of Japanese eel (Chan et al., 1997), the selection pattern between Japanese eel and the North Atlantic eels could be different. Different migration patterns causes distinct selective pressures. With the advances of next-generation sequencing technology, genome wide SNPs have been adopted in population genomic studies of marine fishes. Pujolar et al. (2014) detected signatures of selection using 50 354 RAD-generated SNPs and found dozens of loci associated with the highly variable environmental conditions experienced by European eel along its geographic range. However, study of local adaptation in A. japonica by using genome wide genetic variations is still absent. A better understanding for genetic mechanism of local adaptation in A. japonica could be achieved when genome wide markers were applied in the near future.

5 CONCLUSION

For the first time, our study developed twenty-four gene-associated microsatellite markers and confirmed panmixia of A. japonica by using these microsatellites. The results suggested that A.japonica should be managed as a single unit. Management and conservation efforts must be coordinated at the transnational level, as over-exploitation in any region will decrease the recruitment of this species. Efforts were also made to detect signature of selection and no footprints of spatially varying selection were found, indicating that A. japonica might occupy a different pattern of selection from North Atlantic eels. The study on population genetic structure and local adaptation in A. japonica is far from complete. A whole-genome level study using population genomic approaches is needed in the future, which could help to elucidate the population structure and genetic mechanism underlying the spatially varying selection for geographic populations of A. japonica.

6 DATA AVAILABILITY STATEMENT

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

7 ACKNOWLEDGMENT

We thank Dr. LI Yulong for his help in data analysis. Thanks Dr. LIU Bingjian, Dr. Yokogawa Koji, and Dr. HAN Yusan from National Taiwan University, China, for their help in collecting the samples.