The detailed population genetic structure of the rare endangered latid fish akame Lates japonicus with extremely low genetic diversity revealed from single-nucleotide polymorphisms

Accurate knowledge of spatial population structure is crucial for species conservation and management, but low genetic diversity of endangered species makes it difficult to ascertain this information. In the Japanese lates or akame, Lates japonicus, an endangered fish known to have extremely low genetic diversity among vertebrates, mtDNA and AFLP marker analyses failed to clearly partition individuals into the two main distribution areas (Kochi and Miyazaki prefectures), although a weak but significant genetic difference between these regions was detected. Here, we conducted individual-based clustering analysis and sibship reconstruction to evaluate the species’ population genetic structure by using 282 yearlings collected from four sites in Kochi and Miyazaki prefectures. We used single nucleotide polymorphism (SNP) markers from double-digest restriction-site-associated DNA sequencing (ddRAD-seq). Although our data revealed extremely low genome-wide nucleotide diversity (π = 0.00006) in akame, the clustering analysis showed a clear genetic difference between the two regions, with a small number of migrants and individuals with admixed ancestry. Sibship analysis estimated that the only full-sib pairs of individuals were from the same prefecture and no pairs were from different prefectures. From these results, we conclude that the akame population of Kochi and Miyazaki prefectures are different management units, and for revealing detailed population genetic structure of endangered species which has low genetic diversity, it is effective to conduct individual-based analysis and sibship reconstruction by using juvenile samples and many SNPs from ddRAD-seq.


Introduction
Assessment of population genetic structure and the establishment of appropriate conservation units are essential for conserving rare, endangered species. If genetically distinct local populations exist within a species, each of these populations is considered as an evolutionarily significant unit (Ryder 1986) or a management unit (Moritz 1994), and these units must be conserved separately to prevent loss of the genetic characteristics of each local population (Fraser and Bernatchez 2001). Therefore, for the appropriate conservation of a target species, it is crucial to determine its population genetic structure.
In recent decades, the development of sequencing technology has made it possible for researchers to use genetic markers such as DNA fingerprinting markers and short tandem repeats, and studies using these markers have yielded a lot of useful information for conservation efforts, including population genetic structure (Frankham 2010). Conservation genetic studies using these markers have been conducted in many endangered species from both terrestrial and aquatic environments, and appropriate conservation strategies, including the definition of management units, have been mentioned in these studies (e.g., for Tokyo bittering Tanakia tanago, Kubota et al. 2010; Przewalski's gazelle Procapra przewalskii, Yang and Jiang 2011;yellow cardinal Gubernatrix cristata, Domínguez et al. 2017; and Cuban rock iguana Cyclura nubila nubila, Shaney et al. 2020). However, most endangered species have relatively low genetic diversity because of their small population size . This means that they possess generally less genetic variation in various DNA markers than common species. In addition, obtaining a sufficient number of samples is often difficult for rare endangered species. Such low genetic diversity and sampling difficulty sometimes make it difficult for researchers to determine a species' population genetic structure in detail (e.g., for Siberian salamander Ranodon sibiricus, Chen et al. 2012; hawksbill sea turtle Eretmochelys imbricata, Natoli et al. 2017; and gharial Gavialis gangeticus, Sharma et al. 2021).
A good representation of an endangered aquatic species for which the population structure remains unclear owing to low genetic diversity is the Japanese lates or akame, Lates japonicus. Akame is a large, predatory, endangered marine fish endemic to Japan (Katayama and Taki 1984;Kinoshita and Iwatsuki 1996;Uchida 2005). This species is known to possess extremely low genetic diversity-comparable to that of endangered freshwater fishes which exhibit some of the lowest levels of genetic diversity among vertebrate species (Takahashi et al. 2015)-despite the fact that marine fishes generally have greater genetic diversity than freshwater fishes (McCusker and Bentzen 2010). Its distribution range is quite small, occurring mainly in the estuaries and coastal areas of Kochi and Miyazaki Prefectures in Japan ( Fig. 1), compared with that of the common congeneric species barramundi, Lates calcarifer, which is distributed throughout the Indo-West Pacific region (Greenwood 1976;Katayama and Taki 1984;Iwatsuki et al. 1993). Akame is listed as an endangered species on the Japanese Red List of Japan's Ministry of the Environment (Japanese Ministry of the Environment 2013) and on the Red List of Miyazaki Prefecture (Miyazaki Prefecture 2015), mainly because of the loss of eelgrass (Zostera japonica) bed habitat, which is essential for the fish's larvae and juveniles (Kinoshita et al. 1988;Kinoshita and Iwatsuki 1996). However, in Kochi Prefecture, akame was delisted from that prefecture's Red List in 2017 (Kochi Prefecture 2017), despite a lack of clear evidence of population size recovery and ignoring the conservation genetic evidence mentioned above. The level of conservation strategy is therefore now inconsistent between the fish's two main distribution areas. Furthermore, its range of spawning habitat or spawning areas, which is strongly related to its population genetic structure, remains unclear.
To obtain fundamental genetic information on akame, we previously conducted an amplified fragment length polymorphism (AFLP) analysis (Takahashi et al. 2015). That study showed a weak, but significant genetic difference between Kochi and Miyazaki prefectures, although no detailed information on population genetic structure, including genetic admixture or migration, was obtained from the individual-based analysis. This is mainly because of the extremely low genetic diversity and an insufficient number of genetic markers and samples. In addition, collecting a sufficient number of adult akame from multiple sites in a short period is also difficult because of its rarity. Therefore, further research using many genetic markers and a large number of samples is needed to clarify the population genetic structure of akame for conservation purposes.
In this study, we aimed to clarify the detailed population genetic structure and management units of rare, endangered akame with low genetic diversity and sampling difficulty by using many SNP markers. Rapid progress of sequencing technology has facilitated the large-scale discovery of SNPs, and these SNPs have been applied to conservation genetics (Kumar et al. 2012). We obtained SNPs from double-digest restriction-site-associated DNA sequencing (ddRAD-seq), which provides genome-wide SNPs easily and more cheaply than with other kinds of RAD-seq methods (Peterson et al. 2012). To overcome the difficulty of collecting adult samples, we caught yearling fish in a single cohort from eelgrass bed nurseries where the yearlings aggregate (Kinoshita et al. 1988;Kinoshita and Iwatsuki 1996) and used them for genetic analysis. We conducted the individual-based clustering analysis and sibship reconstruction within yearling samples to examine the detailed population genetic structure. Estimated full siblings from sibship reconstruction analysis were expected to consist of individuals from the same age cohort, meaning that they were born in the same spawning area at the same time. Therefore, we considered that sibship reconstruction might provide direct information about migration distance and spawning area range, thus helping us to understand the population genetic structure. We expected that our results would provide not only the information needed to conserve akame but also useful information for assessing conservation strategies for endangered species with very low genetic diversity and sampling difficulty.  Table 1). All individuals were collected by using a coarse-meshed beach seine (4-mm × 4-mm-square mesh, 10 m wide and 2 m deep with a 3-m-long purse-bag at its center). Histograms of the standard length at each sampling site are shown in Supplementary Fig. 1. Fin clippings were taken noninvasively from all 284 yearling fish, and after the fin clipping and the measurement of standard length all fish were released. This  (Kagimoto and Yamagata 1997;Arai 2005) are shown 50 5 0.2 command. The genotype data set was exported to genepop and structure format by using PGDSpider v 2.1.1.5 (Lischer and Excoffier 2012) and used in further analyses, except in the case of the genetic diversity.

Genetic diversity
We calculated genetic diversity indices at the sampling-site level with all polymorphic loci and all nucleotide positions. To determine the genetic diversity indices, we used loci that were present in at least 80% of individuals and at no fewer than three sampling sites, along with a minor allele count of > 1 in all individuals. Genetic diversity indices, namely the ratio of polymorphic sites, the number of private alleles (A p ), nucleotide diversity (π), expected heterozygosity (H e ), and observed heterozygosity (H o ), were calculated by using Stacks software.

Sibship reconstruction
We used Colony v2.0.6.5 (Jones and Wang 2010) to reconstruct the sibships of sampled individuals in order to evaluate the family bias in juvenile sampling and to obtain information about migration distance and spawning area, which are closely related to the population genetic structure of akame. Colony software analyses were run by assuming polygamy for both sexes, using the full-likelihood method with no sibship prior and with the genotyping error rate set to 0.05. The run length was medium, and the precision of the assignment was very high. Only full siblings that had a probability of assignment above 0.99 were considered as true full siblings. Also, we used the package sequoia (Huisman 2017), implemented in R (R Core Team 2020), and we compared the results of the two methods to validate the accuracy of the estimated sibships. The error rate of the sequoia analysis was set to 0.05, and the parameter MaxSibIter = 10.

Population genetic structure
To evaluate the population genetic structure we used Structure v2.3.4 (Pritchard et al. 2000) The dataset was analyzed by using K (number of clusters) from 1 to 8, with 10 independent iterations, 100,000 Markov chain Monte Carlo (MCMC) repetitions, and a burn-in of 100,000. The most likely number of populations (K) was estimated by using the ΔK method of Evanno et al. (2005) in Structure Harvester (Earl and von Holdt 2012); we determined the optimal K representing the highest ΔK. The software CLUMPP v1.1.2 (Jakobsson and Rosenberg 2007) was used to average the results of the 10 iterations with the identified optimal number of clusters; the results were displayed graphically in order of individual ID in each site by using distruct v1.1 research was conducted with the permission of the Miyazaki and Kochi prefectural governments.

DNA extraction and ddRAD-seq
We used all 192 individuals from K-KMA and M-SHI for genetic analysis. Total genomic DNA was extracted from the fin clippings by using a DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). The concentration of extracted DNA was quantified by using a Quantus fluorometer (Promega, Madison, WI, USA) and a Qubit dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA) and adjusted to 10 ng/ µL. ddRAD-seq was conducted in accordance with the protocol of Peterson et al. (2012), with some modifications. The restriction enzymes EcoRI (6-bp recognition site) and MseI (4-bp recognition site) are frequently used for ddRAD-seq, but we used the restriction enzymes BglII (6-bp recognition site) and MseI to increase the total number of SNPs. Sequencing of 51-bp single-end reads was done by using a Hiseq 2500 sequencing system (Illumina, San Diego, CA, USA) on the 192 samples. Additionally, after 51 bp sequencing, we prepared ddRAD-seq library in the same procedure on all 92 samples from K-SMT and M-OYD to investigate the population genetic structure between sampling sites within the same prefecture. Sequencing of 151bp paired-end reads was done by using a Hiseq X system (Illumina) on these 92 samples.

Genotyping of sequence data
The RAD libraries were demultiplexed by using process_ radtags in Stacks-2.54 (Catchen et al. 2011). Sequence reads were preprocessed by using the Trimmomatic 0.39 read-trimming tool (Bolger et al. 2014), with the following parameters: ILLUMINACLIP TruSeq3-PE-2.fa:2:30:10, LEADING:19, TRAILING:19, SLIDINGWINDOW:30:20, AVGQUAL:20, and MINLEN:51 or 151 (depending on the sequencing length). After the preprocessing, samples with fewer than 1,500,000 reads were excluded (one sample from M-SHI and one from K-KMA), and the remaining reads were mapped to the draft akame reference genome (Hashiguchi et al. in preparation) by using Nextgenmap 0.5.5 (Sedlazeck et al. 2013) with the default settings. Sorted individual BAM files were created with Samtools 1.9 (Li et al. 2009), and SNP calling was achieved with Stacks by using the reference mapping pipeline ref_map.pl in Stacks (with the default parameters). We used PLINK v1.90 (Purcell et al. 2007) to filter out SNPs with minor allele frequency (MAF) < 0.1 for all individuals, missing data > 0.01, and deviation from Hardy-Weinberg equilibrium (HWE) at P < 0.001. Finally, filtered SNPs were pruned for linkage disequilibrium (LD) in PLINK by using the indep-pairwise to archive the recommended acceptance rates between 20% and 60%, as recommended by Wilson and Rannala (2003). We performed 10,000,000 MCMC iterations, with 1,000,000 iterations to discard as burn-in and sampling every 1,000 iterations.
We conducted gene flow model analysis by Fastsimcoal2 software (Excoffier et al. 2013) Supplementary Fig. 9) and 3 datasets from 10 individuals. First, we used the top 5 individuals in terms of sequencing depth from K-SMT, and the top 5 from M-OYD (151-bp paired-end data) to maximize the number of SNPs obtained (dataset A, 3,819 SNPs). We additionally used 2 datasets to evaluate the robustness of model selection and the estimated parameter values. Each top 10 individuals of sequencing depth from K-SMT and M-OYD were divided into two groups (containing 5 individuals from each sampling site) equally in terms of the total amount of sequencing data (dataset B, 3,547 SNPs; dataset C, 3,535 SNPs). Gene flow model analysis based on 100 replicated times with 100,000 simulations for the calculation of the composite likelihood and 40 ECM cycles.

Estimation of population size
The effective number of breeders (N b ) was estimated for the Kochi and Miyazaki populations by using the sibship assignment (SA) method available in the program Colony, as well as the LD method (Waples 2006) by using NeEstimator v2.1 (Do et al. 2014). N b corresponds to the effective population size (N e ) estimated from a single age cohort; it is equivalent to the effective number of adults that contributed offspring to the cohort (Wang 2016). The SA method infers N b from the sibship frequencies with multilocus genotypes of a sample captured at random from a single cohort (Wang 2009). The parameters of estimation were the same as those used in the sibship analysis.

SNP calling and genotyping
The 51-bp single-end sequencing produced a total of 939,361,179 raw reads with an average of 4,892,506 reads across 192 individuals, and the 151-bp paired-end sequencing produced 747,371,582 raw reads, with an average of 8,123,604 reads across 92 individuals. After preprocessing, the remaining reads were mapped to the akame draft (Rosenberg 2004). Other investigations of population genetic structure, namely principal component analysis (PCA) and discrimination analysis of principal components (DAPC), were performed by using the dudi.pca and dapc functions, respectively, in the package adegenet v2.1.1 (Jombart 2008), implemented in R. The find.clusters function was used to determine the number of groups (K), with the optimal K selected as that with the lowest BIC (Bayesian information criterion) value. Cross-validation by using the xvalDapc function in the adegenet package was used to identify the optimal number of principal components to retain in DAPC, with default parameters except for n.rep = 100. We performed a hierarchical analysis of molecular variance (AMOVA) using software Arlequin 3.5.2.2 (Excoffier and Lischer 2010) with 10,000 permutations to reveal the distribution of genetic variation. We calculated the pairwise F ST (fixation index) between each site pair by using Arlequin with 10,000 permutations to detect genetic differences between the sampling sites. The Bonferroni method (Rice 1989) was used to adjust the significance levels when multiple simultaneous tests were performed. Global F ST value was also calculated and tested the significance with 10,000 bootstraps using the package hierfstat (Goudet 2005) in R. Anderson and Dunham (2008) reported that analysis by using Structure and a sample that includes family groups may produce a false population structure, even when a population structure is absent. Therefore, all of the above individual-based population genetic structure analyses were conducted once again after the exclusion of full-sib individuals from Colony to eliminate the effects of family bias due to juvenile sampling on the genetic analyses.

Detection of migrants and recent gene flow
Structure software was also used to identify migrants and those individuals with migrant or mixed ancestry. Sampling location information was used in the analysis to help with the clustering process. We categorized the samples into two groups (Kochi and Miyazaki prefectures) in accordance with their sampling sites. Structure was run with K = 2 from a clustering analysis without population information. As we had no information about migration, the migration rate (MIGPRIOR) was assigned as an initial condition in accordance with the method of Pritchard and Wen (2004). Migrant ancestry (GENSBACK) was set to 1, and the other parameters were the same as for the runs of Structure analysis without population information.
Recent gene flow between the Kochi and Miyazaki populations was assessed by using BayesAss v3.0.4 software (Wilson and Rannala 2003). The MCMC mixing parameter values for migration rate, allele frequency, and inbreeding coefficient were set to 0.08, 0.19, and 0.012, respectively, sites and H e was 0.00006 for all four sites. π was 0.00006 for all four sites (Table 1).

Sibship assignment
Full siblings were identified by using Colony and sequoia. From among the 282 genotyped individuals, Colony assigned 26 individuals into 13 pairs of full siblings that had probabilities of assignment above 0.99. Sequoia assigned 14 individuals into 7 full-sib pairs; these were all included in the Colony result and all had probabilities of 1 (Table 2). Full-sib members were caught at the same sampling sites, or sometimes at different sites in the case of Kochi Prefecture; no full-sib pairs included individuals caught from different prefectures.

Population genetic structure
In the Structure analysis, the results of Structure Harvester showed the highest ΔK at K = 2 ( Supplementary Fig. 2). At genome. Two individuals were excluded because of a deficiency of raw reads, so 282 samples were applied to the Stacks ref_map pipeline; 2,296,333 loci and 92,000 SNPs were initially identified. After filtering for MAF, missing data, HWE, and LD, a total of 780 SNPs were obtained and used for further genetic analysis, with the exception of the genetic diversity and gene flow model analyses. No major batch effects that affect for genetic analysis derived from sequencing length were not found.

Genetic diversity
The proportion of polymorphic sites ranged from 0.0002092  -Inclusion probability: the probability that all individuals of a given full-sib family are full siblings. Exclusion probability: the probability that all individuals of a given full-sib family are full siblings and that no other individuals are individuals within this full-sib family. These probabilities were estimated by using Colony. Assigned by sequoia: ○ indicates that this full sibling was estimated by using sequoia. All full-sib pairs estimated by using sequoia were included in Colony estimations. SL, standard length the population genetic structure analyses after the exclusion of the 26 full siblings assigned by Colony were almost the same as the results obtained by using all individuals (Supplementary Fig. 6 to 8). AMOVA revealed 1.78% of the total variation arose among two prefectures (F CT =0.0178, P < 0.001), 0.23% arose among sampling sites within prefecture (F SC =0.0023, P < 0.001) and remaining 97.99% arose within sampling sites (F ST =0.0201, P < 0.001).
The pairwise F ST values were small between two sites in the same prefecture, at 0.0020 (between K-KMA and K-SMT) and 0.0026 (M-SHI and M-OYD). Between two sites from different prefectures, the pairwise F ST values ranged from 0.0173 (M-SHI and K-KMA) to 0.0240 (M-OYD and K-KMA; M-OYD and K-SMT) (Supplementary Table 2). Global F ST value was 0.0148. All pairwise F ST and global F ST values were significantly different from zero at P < 0.001. K = 2, a distinct genetic difference between the Kochi population and the Miyazaki population was observed, except in the case of some individuals (Fig. 2). At K = 4, no clear difference was apparent among the four sampling sites (Supplementary Fig. 3). We constructed a scatterplot of a PCA analysis by using the first two principal components (PCs) (Fig. 3). The first PC explained 2.9% of the total variation and clearly distinguished the two clusters shown in the Structure analysis. The second PC accounted for 1.1% of the variation. BIC analysis with the find.clusters function revealed K = 2 to be the most likely number of genetic clusters ( Supplementary Fig. 4), and the optimal number of PCs retained for analysis was 20 ( Supplementary Fig. 5). At K = 2, DAPC showed a clear genetic difference between the Kochi and Miyazaki populations, confirming the findings of the Structure and PCA (Fig. 4). The results of all  Plots represent samples and are colored according to sampling site. The x-axis represents PC1 (explains 2.9% of the total genetic variance) and the y-axis represents PC2 (explains 1.1% of the total genetic variance) Fig. 2 Result of a Structure analysis for K = 2, chosen by the Evanno method (Evanno et al. 2005). Each bar shows a different individual, each sampling site is separated by a black line, and the colors represent the probability membership coefficient of each individual for each genetic cluster. ↓ above the bar shows that the individual has admixed ancestry, and ↑ below the bar shows that the individual is a migrant, as estimated from the Structure analysis gene flow) were almost the same with the value of Model b (ΔAIC between Model b and c = 1.326; ΔAIC between Model b and d = 0.074). On the other hand, the best model for dataset B and C was Model d (Supplementary Table 3). The estimated parameters of dataset B and C were quite different although the total SNP numbers of the two datasets were almost the same (Supplementary Table 4). Therefore, we concluded that the results of gene flow model analysis were not reliable.

Estimation of population size
The effective number of breeders (N b ), as determined by using the SA method, was estimated at 459 (95% CI: 363 to 601) in the Kochi population and 376 (95% CI: 301 to 489) in the Miyazaki population. The estimated N b as determined by using the LD method was 551.7 (parametric 95% CI: 522.3 to 597.6) in Kochi and 474.9 (parametric 95% CI: 456.9 to 513.9) in Miyazaki.

Population genetic structure and management units
Our study indicated that genetic analysis using many SNPs enabled us to reveal the detailed population genetic structure of an endangered species that has extremely low genetic diversity. The genome-wide nucleotide diversity (π) of akame was 6.0 × 10 − 5 for all four sites. In similar studies of endangered fishes by using genome-wide SNP markers generated by RAD-seq, the nucleotide diversity of European eel (Anguilla anguilla) was 5.3 × 10 − 3 (Pujolar et al. 2013), that of the salmonid fish Siberian taimen (Hucho taimen) was 1.6 × 10 − 3 (Galland et al. 2021), that of Japanese eel (Anguilla japonica) was 8.0 × 10 − 4 (Gong et al. 2019), and that of Chinese sturgeon (Acipenser sinensis) was 4.8 × 10 − 4 (Liu et al. 2018). The nucleotide diversities of these fishes and some other threatened species of mammals and birds are shown in Table 3. Compared with these values, the nucleotide diversity of akame was extremely low. Although our result might have been underestimated to some extent because of SNP filtering or the selection of conservative regions as a result of application of the ddRAD-seq method, we confirmed the extremely low genetic diversity of akame with SNPs at a genome-wide level. The results of the Structure, PCA, and DAPC analyses using SNPs all showed a distinct genetic difference between the Kochi and Miyazaki populations (Figs. 2, 3 and 4). AMOVA also showed small genetic variation among sampling sites within prefectures (0.23%, F SC =0.0023) compared to the values of among

Detection of migrants and recent gene flow
Structure analysis using population information detected 9 individuals in M-SHI that were migrants (probability of migrants > 0.9), from Kochi to Miyazaki populations (P = 0.987 to 1.000), and 5 individuals as having admixed ancestry, with a migrant parent (probability of having migrant ancestry in one previous generation > 0.9); one of these was caught in Miyazaki (P = 0.997) and four were from Kochi (P = 0.990 to 0.998) (Fig. 2). These 9 migrants had Q-values of > 0.8 (0.8932 to 0.9927) for membership of a non-home cluster in the analysis without population information. Five individuals with admixed ancestry had Q-values in the Structure analysis between 0.2 and 0.8 (0.3371 to 0.5468) for membership of a home cluster, matching the threshold for potential admixture (Lecis et al. 2006;Bergl and Vigilant 2007;Reddy et al. 2012;Ferreira da Silva et al. 2014). Two migrants were assigned as full siblings with each other by Colony.
Recent migration rates between the two populations, as calculated in BayesAss, were 0.0330 from Kochi (K-KMA and K-SMT) to Miyazaki (M-SHI and M-OYD) population, 0.9670 from Miyazaki to Miyazaki, 0.0448 from Miyazaki to Kochi, and 0.9552 from Kochi to Kochi.
Gene flow model analysis revealed that the best model of dataset A was Model b (gene flow only in the direction from Miyazaki to Kochi forward in time), although the AIC value of Model c (gene flow only in the direction from Kochi to Miyazaki forward in time) and Model d (bidirectional Table 3 Levels of nucleotide diversity (π) of akame and other threatened fishes, mammals, and birds calculated from genome-wide SNP markers generated by RAD-seq Common name conducted a Structure analysis. The results of the Structure analysis were almost the same as those of the ddRAD-seq, and the same individuals were estimated to be migrants (data not shown). Individuals that were close in standard length had close individual IDs, because we assigned the IDs roughly in descending order of standard length. It was therefore possible that this could have affected the continuous appearance of migrants in the Structure analysis. The presence of potential migrants and individuals with admixed ancestry indicated that there was ongoing gene flow between the genetically differentiated Kochi and Miyazaki populations, although there is a possibility that some migrants might return to their natal area for spawning. The recent migration rates calculated in BayesAss showed a small but bidirectional gene flow between the two populations, although the level of differentiation of the two populations was too small to reliably estimate migration rates by using this method (Faubet et al. 2007). These results indicate that there are some migrants or individuals with admixed ancestry and that there is therefore a small degree of gene flow between the two populations, but it is not enough to eliminate the genetic difference between them. It is frequently referred that a migration rate of < 0.1 is needed to maintain the demographic dependence of two populations (Hastings 1993;Waples and Gaggiotti 2006;Palsbøll et al. 2006). Indeed, our BayesAss analysis showed that the recent migration rates between the two regions were 0.0330 from Kochi to Miyazaki population, and 0.0448 from Miyazaki to Kochi, this is consistent with the criteria for migration rate. Further analyses using individuals collected from other cohorts will be needed to obtain more information about migrants and individuals with admixed ancestry.
The fact that individuals with admixed ancestry were found in both regions shows that migrants exist in both regions and the direction of transport might be bidirectional. This is consistent with the recent migration rates calculated in BayesAss. The Kuroshio Current flows from southwest to east (the direction from Miyazaki to Kochi) along the Pacific coast of southwestern Japan. A current that splits off from the Kuroshio Current and flows counterclockwise is often observed in the Hyuga-nada Sea, which separates Miyazaki and western Kochi Prefectures (Kagimoto and Yamagata 1997;Arai 2005; Fig. 1). Not only is it possible that there is a small number of adults that actively migrate; these currents may cause passive bidirectional gene flow by transporting akame eggs, larvae, or juveniles between the two regions.

Effectiveness of juvenile sampling
Because of the rarity of adult akame, it is difficult to collect sufficient numbers of adult samples at multiple sites.
prefectures (1.78%, F CT =0.0178). In addition, the assigned full-sib pairs included only individuals caught from a single prefecture; this was consistent with the results of the population genetic structure analysis. From our results, we conclude that akame populations of Kochi and Miyazaki Prefectures are different management units and should be conserved separately, although this would not support the delisting of this species from the red list only in Kochi Prefecture as discussed below. Within each prefecture, no genetic differences between sampling sites were found from STRUCTURE, DAPC or PCA results, AMOVA showed small but significant genetic difference among sampling sites. Further analysis using much more SNPs than this study will verify the existence or nonexistence of fine-scale population structure in the prefectures.
On the other hand, our pairwise F ST values between the two sampling sites from different prefectures were small (F ST =0.0173 to 0.0240) compared to the value from AFLP (F ST =0.06453, Takahashi et al. 2015) although these values were estimated from different genetic markers. They mentioned that the low population divergence in akame might be a result of recent population fragmentation or some degree of gene flow between the two regions. The low level of population divergence in SNPs as well as AFLP and the existence of gene flow from BayesAss and sibship analysis (see below) support that genetic divergence of akame is maintained at a low level due to the gene flow between the two regions.

The existence of migrants and individuals with admixed ancestry
The previous study using AFLP markers failed to detect distinct genetic clusters from individual-based analysis (Takahashi et al. 2015). This study using many SNPs enabled us to reveal the distinct genetic difference between the two regions and the existence of migrants and individuals with admixed ancestry from individual-based clustering analysis, although small genetic difference between the two regions. Structure analyses detected nine individuals as migrants. All of them were caught in M-SHI (in Miyazaki) and none in Kochi, and one individual from Miyazaki and four from Kochi were assigned as individuals with admixed ancestry that had a migrant as a parent. Because most of the migrants were close to each other in terms of individual ID, we suspected that there could have been an artificial influence during the laboratory experiment. For this reason, we conducted a further analysis using these migrants and individuals which has close to migrants in individual ID. Total genomic DNA was newly extracted; we then obtained SNPs by using the genotyping by random amplicon sequencing, Direct (GRAS-Di) method (Hosoya et al. 2019) and has still not been identified, and further research is needed for us to understand the fish's life cycle and preserve its spawning grounds.

Effective number of breeders and effective population size
Effective population size is one of the most important parameters in conservation, because it influences genetic drift, inbreeding, and loss of genetic variation (Frankham 1995;Charlesworth 2009;Franklin 1980) reported that N e = 50 is required to avoid the effects of inbreeding and N e = 500 is necessary to avoid extinction and cope with a changing environment (50/500 rule). Moreover, Frankham et al. (2014) argued that N e > 500 is too low for the retention of evolutionary potential and that a better approximation is N e > 1000.
The relationship of N e and N b has been researched in many species, and in most cases, the relationship is more complicated in iteroparous species than semelparous ones (Waples 2010). In the semelparous fish Chinook salmon, Oncorhynchus tshawytscha, Waples (1990) reported that Ne is equal to gN b , where g is the average generation time. On the other hand, Wang (2009) showed a general relationship of N b ≤ N e ≤ gN b . The N b of the iteroparous fish akame, as estimated here, was 474.9 (LD method) and 376 (SA method) in Miyazaki Prefecture, and 551.7 (LD) and 459 (SA) in Kochi Prefecture. From these estimations, we surmised that the effective population size of akame is 500 to 2000 for each population based on the relationship of N e and N b , taking the average generation time of akame as 5 years. This estimated N e of akame populations of each population may not be sufficient for their persistence for a long time in light of the desired value of N e > 1000, so we will need to conserve the akame populations of Kochi as well as of Miyazaki in appropriate ways. Our results did not exclude the possibility that half-sib pairs were erroneously assigned as full-sib pairs owing, in part, to the extremely low genetic diversity, which leads to the underestimation of N b (by the SA method). The differences in the N b values estimated by using the two methods might reflect such assignment errors. Future studies using SNPs from whole-genome resequencing will help to evaluate sibship assignments and the estimated value of N b of akame with extremely low genetic diversity.

Supplementary Information
The online version contains supplementary material available at https://doi.org/10.1007/s10592-023-01517-2. Therefore, we used many juveniles instead of adults-especially yearling fish, which were collected from eelgrass beds. The standard length distribution of the 284 yearlings collected in this study was 24 to 108 mm, consistent with the criterion about the standard length of yearlings by Uchida (2005). However, yearling samples may contain more closely related individuals than adult samples because of an insufficient dispersal capability of yearlings, and this may lead to bias in allele frequency estimates, which are strongly connected to population genetic structure. In fact, the salmonid brown trout, Salmo trutta, which spawns in limited areas within rivers, has shown significant differences in allele frequencies between age class samples because of family bias in yearling samples (Hansen et al. 1997). The effect of kinship sampling has also been reported in the school shark Galeorhinus galeus (Devloo-Delva et al. 2019). That study showed a significant genetic difference between two groups with the use of neonate or juvenile samples including full-sib groups, whereas in samples without full siblings there was no genetic difference. Our study results without the full-sib individuals were almost the same as those for all individuals including full siblings, showing a distinct genetic difference between the Kochi and Miyazaki populations ( Supplementary Fig. 6 to 8). This indicates that there is no bias for estimating of population genetic structure when using juvenile samples.

Spawning areas of akame
Our research also provides a clue to identifying the spawning areas of akame. The location of these areas has remained unclear for a long time, but there are two main hypotheses. Iwatsuki et al. (1994) and Uchida (2005), on the basis of the reported capture of a female close to spawning at a river mouth, speculated that akame spawn in estuaries or coastal waters. In contrast, Kinoshita et al. (1988), on the basis of a lack of records of mature females in river mouths, inferred that akame do not spawn in estuaries but migrate to the sea and spawn somewhere there. Mitsunaga et al. (2021) proposed that akame in Kochi and Miyazaki Prefectures migrate upstream areas of the Kuroshio Current and spawn around Tanegashima Island, in the Satsunan area of southern Kyushu (Fig. 1). Gonzalvo et al. (2021) conducted an otolith Sr:Ca analysis on juvenile fish caught from the Shimanto River. From their finding of low Sr:Ca ratios in the egg to larval period, they confirmed that akame in the Shimanto River are likely to spawn in the estuary. Our population genetic structure and sibship analyses indicated that each akame population has different spawning areas, thus supporting the estuarine spawning hypothesis and the results of the otolith analysis, rather than the spawning migration hypothesis. However, the exact spawning area