High migratory propensity constitutes a single stock of an exploited cutlassfish species in the Northwest Pacific: A microsatellite approach

Abstract

Cutlassfishes, also known as hairtails, include multiple predatory fishes of the family Trichiuridae. They constitute a top marine fish commodity globally, yet the knowledge about their composition and intraspecific genetic structures is still limited. Trichiurus japonicus accounts for a major amount in the northwest Pacific fishery. Previous studies based on mitochondrial DNA markers reported incongruences in its population structure, hence prompting the need for high-resolution markers and avoiding possible shortcomings in its management. Here we genotyped ten novel de novo-assembled transcriptome-derived microsatellite markers on a total of 150 samples across five major fishing grounds (encompassing latitudes 22–39°N). These markers presented a high number of alleles and heterozygosity compared to other marine fishes, corresponding to the large effective population size of ~20,000 per location and cohort differentiation. Population structuring analyses suggested T. japonicus to be a homogenous well-mixed population. This configuration is likely attributed to the majority of its effective population migrates across locations, and the absence of oceanographic barriers at the continental shelves. Qingdao with reportedly high ocean productivity could be a genetic pseudosink based on the high heterozygosity and migratory preference. Moreover, the results of sign tests suggest that T. japonicus experienced a recent bottleneck likely concurrent with historical glaciation events. Further, we demonstrated satisfactory cross-amplifications of our markers on several congeners, indicating a great promise to use these markers to study the population genetics of trichiurids. Together, our findings will serve as an essential groundwork for enhancing resource conservation and management of cutlassfishes.

About the Authors:

Yu-Hong Guo

Roles Conceptualization, Data curation, Formal analysis, Methodology, Software, Writing – original draft, Writing – review & editing

Affiliation: Department of Marine Biotechnology and Resources, National Sun Yat-sen University, Kaohsiung, Taiwan

Lorenzo C. Halasan

Roles Investigation, Writing – original draft, Writing – review & editing

Affiliation: Department of Marine Biotechnology and Resources, National Sun Yat-sen University, Kaohsiung, Taiwan

ORCID logo https://orcid.org/0000-0002-4752-2231

Hui-Yu Wang

Roles Conceptualization, Data curation, Supervision, Writing – review & editing

Affiliation: Institute of Oceanography, National Taiwan University, Taipei, Taiwan

Hsiu-Chin Lin

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliations Department of Marine Biotechnology and Resources, National Sun Yat-sen University, Kaohsiung, Taiwan, Doctoral Degree Program in Marine Biotechnology, National Sun Yat-sen University, Kaohsiung, Taiwan

ORCID logo https://orcid.org/0000-0001-8116-0484

Introduction

Cutlassfishes or hairtails are predatory fish species of the family Trichiuridae. They are characterized by extremely elongated and compressed bodies and inhabit the tropical and temperate continental shelves [1, 2]. Trichiurids possess high commercial value, which is illustrated by their annual >1.2 million tons catch volume [3]. About 85% of its capture fishery originated from the Northwest (NW) Pacific. In this region, the dominant and most harvested cutlassfish species is Trichiurus japonicus [4–7].

Intra-specific variations were observed in the NW Pacific-endemic T. japonicus. Life-history traits such as growth rate, maturation time, and maximum body length, as well as an inverse correlation to temperature, were found variable [8]. Besides its year-round spawning ability, two spawning peaks observed from the Taiwan coast in February-July and November-December demonstrated the capability of T. japonicus for multiple cohort productions per year [9]. In addition, movements across spawning and feeding grounds, with a close association with the gonadal development, were also detected from the East China Sea observations [10]. These variations were somehow expected for a marine species with a wide latitudinal range (22–39°N). For similar reasons, this also stimulated interests on whether NW Pacific populations constitute a single genetic unit or an integration of genetically different clusters.

Environmental factors such as oceanographic features and geological barriers are known to differentiate populations in various marine organisms [11]. In the NW Pacific, several ocean currents with unique features (e.g. temperature, salinity, and directional flow) interplayed with each other and configured heterogeneous environments, both spatially and temporally. This interplay of current systems is mainly influenced by the Yellow Sea Warm Current (YSWC) [12], the Yellow Sea Cold Water Mass (YSCWM) [13, 14], Taiwan Strait Water Mass (TSW) [15], East China Sea Shelf Surface Water (ECSSSW), and the Kuroshio Branch Currents (KBC) [16]. Furthermore, the massive freshwater discharge of the Yangtze River, also known as the Changjiang Diluted Water (CDW) [17], is also recognized as an oceanographic barrier for some nearby marine organisms [18, 19]. In Tzeng et al. (2016), the partitioning between East and South China Seas populations was attributed to the hindered gene flow of T. japonicus populations, which was caused by the exposure of Taiwan Strait during the glacial periods [20].

Genetic tools have proven to be efficient in identifying and characterizing the presence of different populations or the so-called stocks in fisheries [21–23]. These populations react to natural and anthropogenic impacts differently and thus should be treated as separate units for conservation and sustainable management [24]. Population genetics of T. japonicus have only been studied using mitochondrial DNA markers, including cytochrome b (Cyt-b) [20, 25], control region (D-loop) [26], and 16S rRNA [6] gene regions. Although there is a general pattern of homogeneity and lack of intra-specific differentiation, only one study based on Cyt-b and control region revealed two genetically differentiated populations at the East and South China Seas, with Taiwan Strait as their transient zone [20]. These incongruent results display the limitation of mitochondrial DNA markers in assessing intra-specific patterns for T. japonicus. Other genetic markers, which are non-mitochondrial in origin and have a higher capacity to yield better resolutions, would help elucidate this unresolved intra-specific genetic pattern.

Microsatellites or simple sequence repeats (SSRs) are nuclear DNA with short and tandemly repeated simple sequences. SSRs have been widely used as markers for population genetics due to their high polymorphicity [27]. Their high mutation rate makes them highly variable [28–30], thus providing good genetic resolutions for species targeted for conservation and management [31]. Recent developments in the genotyping protocols and software performances also improved the accuracy of applying microsatellite SSRs [32]. There are two types of microsatellite markers based on where the original sequences were derived from—genomic-SSRs (gSSRs) from genomic sequences and expressed sequence tag-SSRs (EST-SSRs) from transcribed RNA sequences. Ever since the development of next-generation sequencing technologies [33], the recovery of transcribed RNA sequences became cost-effective and less labor-intensive through massively parallel sequencing; and this was then applied to screen for EST-SSRs even for non-model organisms [34, 35]. Compared to gSSRs, EST-SSRs have the advantage of more accurate allele binning [36], high association to functional genes [37], transferability across related species [38], and lesser null alleles [39]. To date, EST-SSRs are still popularly used for studying the population genetics of fishery species [40–42], and some studies have successfully elucidated population structures in the NW Pacific [43, 44]. For cutlassfishes, only a few microsatellite loci have been developed and characterized, and they were exclusively limited to T. haumela [45], T. lepturus [46], T. nanhaiensis [47], Eupleurogrammus muticus [48], and Lepturacanthus savala [49]. Only one study subsequently applied them to conduct population genetic analyses in the Savalai hairtail (L. savala) along the coast of China [50]. Population genetic studies using microsatellite markers on T. japonicus from the NW Pacific are still unavailable.

This is the first study to use microsatellite markers in determining the population genetics of T. japonicus. We hypothesized the presence of an isolated population in the Yellow Sea based on its semi-enclosed topology, unique oceanographic dynamics, and the freshwater barrier from the Yangtze River (Fig 1), as well as in Taiwan Strait where historical vicariance events occurred. We performed transcriptome sequencing using Illumina and developed microsatellite markers accordingly. These markers were then used to study the population structure and demographic history of T. japonicus, which provide relevant information to determine the genetic boundaries for management and formulate sustainable cutlassfish fishery resource strategies.

[Figure omitted. See PDF.]

Fig 1. Map of the Northwest Pacific illustrating the sampling locations of T. japonicus for microsatellite study (circles) and RNA extraction (square).

The details of sampling location and time are described in Table 1. The map was generated by Schlitzer, Reiner, Ocean Data View, odv.awi.de, 2021.

https://doi.org/10.1371/journal.pone.0265548.g001

Materials and methods

DNA sample collection and extraction

Thirty adult T. japonicus individuals were randomly selected from the fisheries catch at each of the five fishing grounds covering the distribution range along the NW Pacific coast, namely: Dalian (DL), Qingdao (QD), Zhoushan (ZH), Badouzi (GE), and Tsukuan (T), from August 2018 to November 2019 (Fig 1, Table 1). A muscle tissue sample from each individual was isolated and preserved in 95% Ethanol. Total DNA used in the microsatellite experiments was extracted from the tissue sample using BioKit Tissue & Cell Genomic DNA Purification Kit (BioKit, Taiwan) following the manufacturer’s protocol.

[Figure omitted. See PDF.]

Table 1. Sampling location and time of T. japonicus specimens collected in this study.

https://doi.org/10.1371/journal.pone.0265548.t001

RNA sample collection, extraction, and sequencing

One T. japonicus individual was collected off Guishan Island, Taiwan (Fig 1, Table 1) in February 2020 for RNA extraction. A fresh muscle tissue sample was preserved in Ambion RNAlater solution (ThermoFisher Scientific, USA) prior to the extraction. Total RNA was extracted with RNeasy Mini Kit (Qiagen, Germany) and eluted with 40 μl of RNase-free water from the RNeasy Mini columns. Following the extraction, the purity of the RNA was assessed and quantified using NanoDrop 2000 spectrophotometer (ThermoFisher Scientific, USA). Fragments of the RNA were also assessed in a 1.5% agarose gel electrophoresis. After the quality control procedures, the mRNA was enriched from total RNA using oligo (dT) beads and fragmented randomly in a fragmentation buffer. Random hexamers were then used to reverse-transcribe the mRNA into cDNA. The cDNA library of T. japonicus sample was prepared by Illumina TruSeq RNA Sample Prep Kits v2 and sequenced using HiSeq 2500 High-Throughput Mode v4 with paired-end 125 basepair reads operated by Novogene Co., Ltd., according to the manufacturer’s instructions (Illumina, San Diego, CA).

Sequence de novo assembly

The sequence read quality was checked and filtered with FastQC v0.11.5 [51] and subjected to Trimmomatic v0.35 [52] for adaptor trimming and rejection of low-quality reads. In obtaining the assembled transcriptome, the clean reads were de novo assembled in Trinity v2.2.0 [53, 54]. Corset v1.06 [55] was then used to remove contig redundancies by hierarchical clustering. Afterward, the longest transcripts of each cluster were selected as unigenes and were used for the succeeding microsatellite loci detection.

Microsatellite locus analysis and primer design

Microsatellite identification tool (MISA) v1.0 was used for microsatellite loci detection in the unigenes [56, 57]. The following minimum number of repeats of each unit size were established as default parameters: mononucleotide repeats ≥ 10; dinucleotide repeats ≥ 6; tri-, tetra-, penta- and hexa-nucleotide repeats ≥ 5. Studies have shown that more information can be drawn out from higher repeat numbers, whereas repeat numbers <10 have already low variability [58, 59]. Perfect repeats were also observed to be more polymorphic than imperfect repeats [60]. Therefore, we did preliminary tests on microsatellite loci with perfect repeats; and the top 22 highest number of repeats, which all have >10 repeats, were utilized as markers for this study. The sequences generated were used to design primers with the Primer3 v2.3.5 program [61] based on the following criteria: (1) primer lengths must be 18 to 23 bp and (2) polymerase chain reaction (PCR) product size must be below 350 bp. Afterward, 22 pairs of primers were manufactured and tested for PCR amplification. Based on initial screenings, 10 pairs ensued successful amplification and thus were chosen as the final microsatellite markers for the succeeding analyses (Table 2). The 5’ end of each forward primer was modified to contain fluorescence (FAM-, HEX-, NED-, and TET-).

[Figure omitted. See PDF.]

Table 2. The characteristics of the ten microsatellite loci of T. japonicus developed from the transcriptome dataset.

https://doi.org/10.1371/journal.pone.0265548.t002

Microsatellite loci genotyping and neutrality test

All T. japonicus samples were genotyped for the ten microsatellite loci: TJ-2, TJ-7, TJ-8, TJ-9, TJ-10, TJ-14, TJ-17, TJ-18, TJ-20, TJ-21 (Table 2). Microsatellite loci were amplified by PCR with primers designed above. Two PCR amplification conditions were applied: (1) For TJ-7, TJ-8, TJ-9, TJ-10, TJ-14, TJ-18, and TJ-21, the solution contained 10 μl Taq 2x Master Mix red (4 mM MgCl2) (Ampliqon, Denmark), 1 μl of template DNA, 0.5 μl of each primer, and ddH2O with a final volume of 20 μl; (2) For TJ-2, TJ-17, and TJ-20, the solution contained 10 μl SuperRed PCR Master Mix (2x) (4 mM MgCl2) (BIOTOOLS Co., Ltd, Taiwan), 1 μl of template DNA, 0.8 μl of each primer, and ddH2O with a final volume of 20 μl. Reactions had an initial step at 95°C for 3 min, followed by 40 cycles of denaturation at 94°C for 30 s, annealing at 54.5–66.5°C depending on the primers (Table 2) for 30 s and extension at 75°C for 45 s, and a final extension at 75°C for 5 min. Microsatellite DNA genotyping was performed using the 3730XL DNA Analyzer (Genomics, Taiwan). Genotypes were confirmed at Peak scanner v2.0 (ThermoFisher Scientific, USA), and the output microsatellite datasets were manually evaluated following Daniels et al. (1998) [62]. These datasets were then tested for genotyping errors, such as stuttering, allelic dropout, and presence of null alleles using MICRO-CHECKER v2.2.3 [63].

The departure of each microsatellite locus from Hardy-Weinberg equilibrium was assessed following procedures described by Louis and Dempster (1987) [64] and performed in GENEPOP [65, 66], with the p-value estimated by Markov Chain (MC) algorithm [67]. The default settings of the parameters were: dememorization number = 1000, batches = 100, and iterations = 1000. Detection of loci under selection based on F-statistics was performed in Arlequin 3.5 [68]. Coalescent simulations were used to obtain the p-values of locus-specific F-statistics conditioned on observed levels of heterozygosities [69, 70]. Unusual genetic differentiation levels will show if loci are detected under natural selection [69].

Genetic diversity

Expected heterozygosity (HE) and observed heterozygosity (HO) were calculated for each sampling location following Weir and Cockerham (1984) [71] using GENETIX v4.05 [72]. Observed heterozygosity (HO) was directly calculated from the proportion of heterozygotes in the sample, whereas expected heterozygosity (HE) was calculated following Nei (1978) [73]. GENETIX v4.05 was also used to calculate fixation indices (FST) and inbreeding coefficients (FIS) between and within each location based on the F-statistics method [71] with 1000 bootstrap replicates.

Population structure

The potential genetic differentiation of T. japonicus was assessed through four methods. First, an unrooted neighbor-joining (NJ) phylogenetic tree was constructed for the five locations using the DA distance method [74] with 1000 bootstrap replicates by POPTREE2 [75]. Second, Bayesian clustering was performed in STRUCTURE v.2.3.4 [76]. A burn-in period of 100,000× followed by 1,000,000× iterations on the Markov Chain Monte Carlo (MCMC) simulations was carried out. The number of separate clusters (K) was set from 1 to 5 and 10× iterations for each specific K-value. Then, the most probable K value based on the posterior probability was determined and subsequently used to assign each T. japonicus individual to its most probable cluster by Structure Harvester [77]. An ad hoc statistic of ΔK values was also calculated based on the rate of change of log probability of the data [78]. Third, Principal Coordinate Analysis (PCoA) was conducted through GenAlEx v.6.503 [79]. The inter-individual genetic distance of microsatellites loci was calculated following Smouse and Peakall (1999) [80] on the methods for codominant data. Fourth, the global and hierarchical partitioning of genetic structures were tested using the Analysis of Molecular Variance (AMOVA), setting codominant genotypic data for distance calculations, by GenAlEx with 9,999 permutations. Based on oceanographic features in the NW Pacific and the results of STRUCTURE v.2.3.4, we classified four types of hypothetical groupings for the five locations in the AMOVA analyses: (1) Five regions (DL; QD; ZH; GE; T); (2) Two regions (Yellow Sea: DL/QD; East China Sea and Taiwan Strait: ZH/GE/T); (3) Three regions (Yellow Sea: DL/QD; East China Sea: ZH/GE; Taiwan Strait: T); (4) Two regions (STRUCTURE results: QD; DL/ZH/GE/T).

Effective population size, gene flow, and bottleneck effect

The effective population size (Ne) and the number of migrants per generation were calculated based on the Maximum likelihood (ML) method in MIGRATE v1.2.32 [81]. Ne was inferred from the frequency of allele exchange among the locations and was calculated by the estimated population size parameter (mutation rate = 10−3) [82]. The assumption of a recently declined population, also known as a ‘bottleneck effect’, states that allele diversity reduces faster than heterozygosity due to the rapid loss of rare alleles [83, 84]. Compared to the expected heterozygosity in a population of constant size with the same number of alleles, a bottleneck effect signifies a temporary excess of heterozygosity [85, 86]. We used two tests to determine which population exhibits a significant number of loci with excess heterozygosity, namely the sign test [85] and Wilcoxon signed-rank test [87], which were both implemented in Bottleneck 1.2.02 [86]. The Wilcoxon signed-rank test is recognized as more powerful than the sign test and can produce better results for fewer polymorphic loci (n = 10 to 15) and datasets with a lower number of individuals (n = 15 to 40) [85]. The simulation result on the population containing a similar number of alleles was compared to the observed heterozygosity of a population evolving under the Infinite Allele Model (IAM) [88] and Two-Phase Model (TPM) [89]. The proportion of the Stepwise Mutation Model (SMM) in TPM was set to 5% and run for 100,000 iterations. The TPM model was hypothesized to fit with the microsatellite DNA data than that of the IAM model due to parameter explicitness [86, 90].

Application of microsatellite markers in congeneric species

The application of the ten microsatellite markers developed in this study was also tested on three congeneric species, T. brevis, T. lepturus, and T. nanhaiensis, which also occur in the NW Pacific waters. Two DNA samples of each species were collected, extracted, and subjected to the PCR with conditions elaborated earlier.

Results

De novo RNA assembly and microsatellite loci identification

A total of 47,923,043 clean reads were generated after adaptor removal from the Illumina sequencing products. The length of assembled sequences ranged from 301 bp to 66,280 bp, with an average of 1,776 bp. GC content was 51.11%, and the length of N50 was 2,852 bp. Ensuing transcriptome assembly, a total of 46,632 assembled unigenes were identified. The microsatellite analysis using MISA examined the unigenes of T. japonicus transcriptome and found 40,027 (85.84%) that contained microsatellites. On all the motif repeat sizes, mononucleotide was the most abundant (41.29%) (S1 Table). A total of 10 pairs of primers designed in this study successfully amplified the targeted microsatellite regions, including 9 pairs for dinucleotide repeats and one pair for trinucleotide repeats (Table 2).

Microsatellite genotyping and neutrality test

Genotyping errors caused by the presence of stutter and allele dropouts were not detected by MICRO-CHECKER. Null alleles were potentially present at the locus of TJ-2 in ZH and GE; TJ-9 in GE and T; TJ-10 in QD and GE; TJ-14 in ZH; TJ-18 in ZH and GE; and TJ-21 in GE and T, as suggested by the general excess of homozygotes for most allele size classes. Almost all loci showed a significant departure from Hardy-Weinberg equilibrium (P < 0.05) (S2 Table). Detection screening of loci under selection from F-statistics showed that none of the ten loci had unusual genetic differentiation levels, thus no sign of selection (S1 Fig). As the result, genotyping data of all ten loci were included for the following analyses (S3 Table).

Genetic diversity

The number of alleles (Na), expected heterozygosity (HE), and observed heterozygosity (HO), per locus and location, were listed in S4 Table. The number of alleles ranged from 10 (TJ-17, ZH) to 36 (TJ-2, ZH and T), with an average of 23.14. Of the ten loci, TJ-2 (33.6) and TJ-17 (13.6) had the highest and lowest number of alleles, respectively. Of the five locations, GE (23.5) and ZH (22.9) had the highest and lowest number of alleles, respectively. The expected heterozygosity ranged from 0.8461 (TJ-17, ZH) to 0.9625 (TJ-2, ZH). The mean expected heterozygosity per locus was the highest in TJ-2 (0.95656) and lowest in TJ-17 (0.86978). The mean expected heterozygosity per location was highest at T (0.9253) while lowest at ZH (0.9156). The observed heterozygosity ranged from 0.7241 (TJ-18, ZH) to 0.9667 (TJ-7, QD). The mean observed heterozygosity per locus was the highest in TJ-7 (0.92), and the lowest in TJ-20 (0.84666). The mean observed heterozygosity per location was highest at QD (0.844) and lowest at GE (0.7939).

The F-statistics results of FST and FIS were shown in S5 Table. The FST values ranged from -0.0043 (TJ-17, QD) to 0.01476 (TJ-18, T). For the FST values of the ten loci across all locations, the highest value was recorded from ZH (0.008019), while the lowest value was from DL (0.002937). All FST values were below 0.05, which indicated that there is no differentiation in each location. The FIS values ranged from 0.00819 (TJ-14, QD) to 0.13604 (TJ-9, GE). For the FIS values of ten loci in all five locations, the highest value was observed at DL (0.072342) and lowest at ZH (0.065282). All the FIS values and the 95% confidence interval (0.04785–0.08527) were very low, indicating a low level of inbreeding in each location.

Population structure

All four methods suggested a mixed population for T. japonicus. The NJ tree showed poor resolution on the genetic relationships of the five populations of T. japonicus, represented with low node supporting values (Fig 2). STRUCTURE Harvester showed that ΔK was highest at K = 2 (Fig 3A), with the most probable K value at K = 1 and 2 (Fig 3B). The result of K = 1 suggested that the five locations were not further separated into different clusters. Meanwhile, the result of K = 2 showed that the majority of individuals from all locations comprise a large single cluster, while few individuals from QD are categorized into a smaller second cluster (Fig 3C). The PCoA result based on co-dominant genotypic distances revealed that all individuals from five locations are mixed (Fig 4). The first three axes accounted for 7.99% of the total variation, and 2.94%, 2.56%, 2.49% for each. The pairwise FST values among locations ranged from 0.00148 (between GE and T) to 0.00766 (between DL and QD). They were significantly different for DL/QD/ZH vs T; QD vs DL; ZH vs T; and ZH vs the rest of the locations (S6 Table). AMOVA results indicated that there was no differentiation among locations under the four types of hypothetical groupings. The percentages of total variation among regions were all less than 1%, and most of the genetic variations were within locations with percentages higher than 98% (S7 Table).

[Figure omitted. See PDF.]

Fig 2. NJ tree based on ten microsatellite dataset from five T. japonicus populations.

The scale bar at the bottom indicates the number of expected substitutions per site. Numbers by branch are bootstrap values.

https://doi.org/10.1371/journal.pone.0265548.g002

[Figure omitted. See PDF.]

Fig 3. Hierarchical STRUCTURE analysis based on ten microsatellite datasets from 150 T. japonicus collected from five locations along the coasts of the NW Pacific.

ΔK (A) and ln P(D) (B) are plotted against the number of genetic clusters (K). The population structure result when K = 2 is shown in (C). Each vertical colored line represents the assignment proportion of membership of each sample and categorized by locations.

https://doi.org/10.1371/journal.pone.0265548.g003

[Figure omitted. See PDF.]

Fig 4. Plot of the first two axes of Principal Coordinate Analysis (PCoA) based on codominant genotypic distances of ten microsatellite loci for 150 T. japonicus individuals.

Each dot represents an individual and is color-coded for the five locations reflecting latitudinal gradient (north: blue, south: red).

https://doi.org/10.1371/journal.pone.0265548.g004

Effective population size, gene flow, and bottleneck effect

The effective population size (Ne) of DL, QD, ZH, GE, and T were 21,343.79; 24,014.94; 19,529.13; 20,971.69; and 22,655.89, respectively. The number of migrants per generation ranged from 15,739.67 (DL to T) to 19,656.04 (GE to QD) (S8 Table)—more than half of the effective population size per generation—which showed significant migration activities between locations. Upon ranking these migrant estimates, QD appeared to be the highly-migrated sampling location by all the other populations. Also, all ten loci fit both the IAM and TPM mutation models. According to the sign test, all locations except QD had significant heterozygosity excess under the IAM model; only T had significant heterozygosity excess under the TPM model. According to the Wilcoxon signed-rank test, all locations had significant heterozygosity excess under the IAM model; all locations except QD had significant heterozygosity excess under the TPM model (Table 3).

[Figure omitted. See PDF.]

Table 3. Bottleneck effect tests for five populations of T. japonicus.

https://doi.org/10.1371/journal.pone.0265548.t003

Application of microsatellite markers in congeneric species

The ten microsatellite markers developed from T. japonicus were highly successful in cross-amplifying with the three congeneric species (S2 Fig). For all three species, TJ-2, -7, -10, -14, -20, and -21 were most promising by showing two clear bands around the desired PCR size range in the gel photo. Additionally, TJ-8, -9, and -18 were promising exclusively for T. lepturus and T. nanhaiensis. However, some adjustments in PCR conditions (i.e., screening for optimum annealing temperature) were required for amplifying the pure target products in some cases.

Discussion

EST-microsatellite markers

Genome-wide genetic information can help us better understand the biological characteristics and molecular mechanisms of organisms. In recent years, Next Generation Sequencing (NGS) has become an efficient and cost-effective method to generate genome-scale sequence data, which advanced the formulation of molecular markers for both model and non-model organisms [91]. Cutlassfish is a non-model organism, and its genome data is currently unavailable. In this study, we generated the first cutlassfish transcriptome dataset from the highly exploited T. japonicus in the NW Pacific using one of the NGS platforms—Illumina, which will lay the foundation for future related studies.

Molecular markers are an important tool to study population genetic diversity. Microsatellite DNA markers are particularly useful in monitoring genetic variability due to their inherent high mutation rates, thus have also been successfully applied to the management of fisheries species [92]. In this study, we identified 40,027 microsatellite markers from the transcriptome dataset. Similar to other fish species, the high proportion of identified microsatellite DNA markers (41.29%) were mononucleotide repeats [40, 93–95].

The null hypothesis on genetic studies states that populations are selectively neutral, meaning they are not influenced by natural selection but by random fixation [96, 97]. Therefore, the departure from Hardy-Weinberg equilibrium of the loci due to selection might not be an ideal marker for population genetic studies. Out of the final ten microsatellite loci applied in this study, some significantly departed from Hardy-Weinberg equilibrium due to heterozygote deficiencies. Probable reasons could include genetic drift, inbreeding, selection, and the presence of null alleles. However, these probabilities were ruled out due to a lack of justification. Since genetic drift often has more impacts on small populations [98], this becomes unlikely for T. japonicus being the dominant cutlassfish species in the NW Pacific and the estimated effective population size per sampling locality also ranged between 19,529 to 24,014. Inbreeding is neither a feasible reason because the inbreeding coefficients are positive yet very low [99]. We did not see evidence of either natural or human-induced selection, such as fishing pressure, based on the F-statistics results. A microsatellite null allele can be derived from a failed amplification during the PCR process [100]. In this study, some potential null alleles were suggested yet their homozygous identity was validated by reperforming the PCR using lowered annealing temperature. Carlsson (2008) also showed that despite the presence of null alleles, this does not necessarily contribute a big impact on the subsequent analyses, compared to other factors such as the strength of population differentiation and the number of loci [100].

Another potential explanation for the departures from Hardy-Weinberg equilibrium is the Wahlund effect [101]. In population genetics, Wahlund effect is the reduction of heterozygosity in a population due to the presence of subpopulation separations by geographic barriers or age structure [102]. Cutlassfish in the NW Pacific have multiple spawning seasons and multiple spawning sites [9, 103] contributing to temporal and spatial subpopulations [102, 104]. If two or more subpopulations have unique allele frequencies, the overall heterozygosity would be lower than expected—assuming there is a random mating population [101]. Wahlund effect should be expected for the mixed-age population at individual loci [102], such as in T. japonicus.

Finally, we tried to remove loci with departures from the Hardy-Weinberg equilibrium from multiple locations, repeated all analyses, and still obtained the same results. Therefore, we consider that these non-neutral loci might have a weak influence on our subsequent analyses. Also, the inclusion of multiple loci would help average the effects of selection [105]; therefore, all the ten microsatellite markers were retained for the succeeding population genetic analyses.

High genetic diversity

Genetic diversity indices, including the number of alleles and observed heterozygosities of the ten microsatellite markers showed higher values in T. japonicus than other reports on commercial marine fishes. The average number of alleles per locus of T. japonicus (23.14) was relatively higher than that of striped bass Morone saxatilis (14.125) [42], pompano Trachinotus ovatus (5.875) [93], and small yellow croaker Larimichthys polyactis (5.82) [106]. Whether these differences might be species or marker-related, the average number of alleles per locus of T. japonicus was relatively higher than 12 other marine fish species (20.6) based on 66 loci [107]. The average observed heterozygosity per locus (0.82514) also showed the same findings, which is higher than other commercial fishes, such as striped bass Morone saxatilis (0.766) [42], pompano Trachinotus ovatus (0.318) [93], and yellow drum Nibea albiflora (0.8033) [108]; and also higher than the average of marine fishes (0.79) [107]. Among the published microsatellite genetic parameters of cutlassfish species, the average number of alleles per locus (23.14) and the average observed heterozygosity per locus (0.82514) of T. japonicus were higher than other cutlassfish species such as T. haumela (Na = 4.58, Ho = 0.661; Na = 12.3, Ho = 0.71) [45, 109], T. lepturus (Na = 11.6, Ho = 0.804) [46], T. nanhaiensis (Na = 4.4, Ho = 0.807) [47], and L. savala (Na = 9.11, Ho = 0.634; Na = 12.5, Ho = 0.81) [49, 50].

The high genetic diversity indices typically indicate a larger population size and/or higher mutation rate [110]. Since mutation rates should be similar among closely related species [111], a larger population size could be responsible for the relatively higher genetic diversity of T. japonicus. In the NW Pacific, T. japonicus has the widest distribution range compared to its congeners [6] and has been reported as the dominant cutlassfish species [7, 112], with an estimated effective population size per location of around ~20,000. In addition, Hsu et al. (2007) proposed that the high genotypic differentiation among cohorts is the underlying mechanism for the genetic diversity in cutlassfish [113]. This differentiation was contributed by the high fecundities and larval mortalities of the species, which in response increase the variance in reproductive success for different cohorts [114]. This heterozygosity is commonly used as an indicator of the informational power of a marker. The heterozygosity values per locus ranged from 0.84666 to 0.92 for the ten microsatellite markers designed in this study, thus indicating high polymorphicity and ideal for population genetic analyses [115].

A well-mixed population experienced bottleneck effect

The results of the NJ tree, AMOVA, PCoA, and STRUCTURE all suggested a single well-mixed population of T. japonicus in the NW Pacific. The NJ tree yielded poor phylogenetic resolution to separate populations, which was indicative of a homogenous structuring. AMOVA showed that there was no significant differentiation at all geographic hierarchical levels. PCoA displayed that all samples were mixed in the quadrant and no obvious grouping was evident. The results of STRUCTURE showed that both K = 1 and K = 2 population schemes were the most probable; however, only a few individuals from Qingdao were assigned to the second cluster, and the reason behind that is still unknown. Although differentiation among some locations was statistically significant based on F-statistics, the FST values were all much lower than 0.05 suggesting minimal differentiation [99]. In addition, the number of migrants per generation among locations was more than half of the effective population size which can greatly facilitate the admixture. The slightly higher migration preference towards Qingdao might explain its relatively high observed heterozygosity. Since populations prefer migrating to Qingdao, this location might have become a genetic pseudosink. Pseudosinks are characterized as viable habitats, particularly due to the very high immigration rate [116, 117]. We suspect that this migratory preference might be related to habitat productivity. Qingdao particularly lies adjacent to the Bohai Sea-Yellow Sea boundary. Nutrient inputs carried by the intrusive YSWC might be partially stalled at the narrow Bohai Sea opening and promotes habitat productivity. In addition, the nutrients discharged in the Bohai Sea and Huanghe River, flushed through by the Bohai Sea Coastal Current also girds around the Shandong coast [118, 119]–where Qingdao is particularly located. This dynamic potentially makes Qingdao an ideally productive foraging ground as reflected also by numerous eutrophication reports [119–121]. This overall high dispersal behavior of T. japonicus was also previously reported to be seasonal in the China Seas and potentially related to spawning seasons [5]. Our finding of a single well-mixed population of T. japonicus is congruent with previous genetic results using mitochondrial DNA markers [25, 26]. This is likely attributed to the absence of geographic isolation in the NW Pacific and the great migratory ability of this fish. Marine fishes can have low levels of genetic differentiation among locations because of the relatively high dispersal potential during pelagic egg, larval, or even adult stages and lack of physical barriers to the movement [122, 123].

A recent bottleneck event was proposed for the majority, if not all locations, based on the sign tests. Although the timing of this event was not estimated, we suspected that it corresponded to the sea-level change in the Pleistocene. It had been well recognized that the cyclical fluctuations of the global climate in the Pleistocene caused drastic changes in sea level, which can affect the distribution range and biomass of marine organisms [123, 124]. This sea-level change can shape the population history of marine species [25, 125, 126]. In addition, biological characteristics of species such as dispersal ability and oceanic characteristics of the environment also contribute to the genetic structure of marine organisms [127, 128]. During the glacial period, the coastline was pushed offshore, and marine residents were forced to new geographic areas or refugia. For continental shelf species in the NW Pacific such as cutlassfishes, suitable habitat was reduced to an elongated enclosed sea in the East China Sea and the South China Sea (the Okinawa Trough) with an area <1/3 of its present size [129], which served as a refuge for the survivors [130].

The Bohai Bay, the Yellow Sea, and the Taiwan Strait all disappeared completely. After the glacial period, rising sea levels brought the fauna back to a broader habitat and open opportunities to expand the population size. He, et al. (2014) and Xiao, et al. (2014) proposed a population expansion of T. japonicus in the late Pleistocene based on mitochondrial DNA markers [25, 26]. The great dispersal ability of cutlassfish and the absence of prominent oceanographic barriers in this area greatly facilitate the expansion process. This also explains the lack of significant population genetic differentiation in the continental shelf of China, which has been reported in taxon-wide marine organisms, including shrimps, polychaete worms, gastropods, and fish [130].

Fishery management

A downward trend of the total catch of cutlassfish in the NW Pacific area has been reported [131–134]. Some studies suggested that persisting fishing pressure had changed their growth and reproduction pattern, as seen in their precocious maturity and miniaturization [5, 8, 135], and this led to fitness reduction [136]. In this study, we found that T. japonicus has high genetic diversity and large population size, suggesting that this highly exploited species is resilient to previous impacts. Moreover, a single population with high gene flow across locations had been identified throughout the NW Pacific. Recent fisheries data have shown that the top exploited Trichiurid relative, T. lepturus, has already been experiencing a decreasing catch rate since the last decade [132]. This serves as a warning and calls for urgent T. japonicus management schemes as this species could potentially be the next sought target.

Based on the above genetic characteristics, we suggest recognizing the entire sampling localities as a single fishery unit and developing international cooperative management. More stringent measures can also be applied to fishing grounds with relatively higher recorded heterozygosity and migration preference (i.e., Qingdao) to establish more substantial protection for the existing genetic diversity. The successful amplification of the EST-microsatellite markers developed in this study towards congeneric species implies their usefulness for characterizing genetic diversity on other highly targeted Trichiurids and confirms the advantage of transferability across related species which economized time and resources [38]. With these cross-amplified markers, it can pave the way for the establishment of genetic marker repositories for the vulnerable Trichiurid species from the NW Pacific.

Conclusions

This is the first study applying microsatellite markers to depict the population genetics of T. japonicus. The ten microsatellite loci developed from transcriptome displayed high polymorphism and genetic diversity thus proved as ideal markers to study intra-specific variations. We found that T. japonicus has higher genetic diversity than other marine fishes including the cutlassfish relatives, which might be attributed to its large population size and/or cohort genetic differentiation. Besides, a well-mixed single population was identified for T. japonicus in the NW Pacific, which has experienced a bottleneck and subsequent expansion likely facilitated by the strong migratory ability and the lack of dispersal barriers. High heterozygosity and migration preference were also reflected for Qingdao, and hence considered as a genetic pseudosink. Based on the above genetic characteristics, we suggest treating the entire NW Pacific population as one fishery unit and developing international cooperative management for the species. The developed microsatellite markers have proven effectiveness for transferability in congeneric species investigations.

Supporting information

S1 Fig. FST (small circles) for individual microsatellite locus plotted against heterozygosity.

Dashed lines are one-sided confidence interval limits obtained from simulation-based expected distributions of FST assuming a stepwise mutation mode.

https://doi.org/10.1371/journal.pone.0265548.s001

(DOCX)

S2 Fig. Gel photo of the PCR results of the ten microsatellite markers in four Trichiurus species (Tj: T. japonicus, Tb: T. brevis, Tl: T. lepturus, and Tn: T. nanhaiensis).

https://doi.org/10.1371/journal.pone.0265548.s002

(DOCX)

S1 Table. The total number and percentage of simple sequence repeats (SSRs) identified for each motif length.

https://doi.org/10.1371/journal.pone.0265548.s003

(DOCX)

S2 Table. The results of Hardy-Weinberg equilibrium analyses for each locus in five locations.

https://doi.org/10.1371/journal.pone.0265548.s004

(DOCX)

S3 Table. The genotyping results of each individual.

https://doi.org/10.1371/journal.pone.0265548.s005

(XLSX)

S4 Table. Genetic diversity estimated based on the ten microsatellite loci in five populations of T. japonicus.

Na: number of alleles; He: expected heterozygosity; Ho: observed heterozygosity. Dark gray means number of alleles > 30 and heterozygosity > 0.9; medium gray means number of alleles between 20 to 30 and heterozygosity between 0.85 to 0.9; light grey means number of alleles < 20 and heterozygosity < 0.85.

https://doi.org/10.1371/journal.pone.0265548.s006

(DOCX)

S5 Table. The F-statistics index estimated based on the ten microsatellite loci in five populations of T. japonicus.

FST: The fixation index; FIS: The Wright’s fixation index.

https://doi.org/10.1371/journal.pone.0265548.s007

(DOCX)

S6 Table. Matrix of pairwise FST among five populations based on T. japonicus microsatellite data.

https://doi.org/10.1371/journal.pone.0265548.s008

(DOCX)

S7 Table. Analysis of molecular variance (AMOVA) results for the four types of groupings.

d.f., degree of freedom.

https://doi.org/10.1371/journal.pone.0265548.s009

(DOCX)

S8 Table. The number of migrants per generation among locations.

Nm = Number of migrants. Location 1 = DL, 2 = QD, 3 = ZH, 4 = GE, and 5 = T.

https://doi.org/10.1371/journal.pone.0265548.s010

(DOCX)

S1 Raw images.

https://doi.org/10.1371/journal.pone.0265548.s011

(PDF)

Acknowledgments

The authors would like to thank Chia-Hong Liu for assistance in sample collection; Jen-Wei Kuo for assistance in conducting molecular experiments and preparing figures and tables; ACJ Hurboda for English grammar checking and additional proofreading.

Citation: Guo Y-H, Halasan LC, Wang H-Y, Lin H-C (2022) High migratory propensity constitutes a single stock of an exploited cutlassfish species in the Northwest Pacific: A microsatellite approach. PLoS ONE 17(3): e0265548. https://doi.org/10.1371/journal.pone.0265548

References

1. Nakamura I, Parin NV. Snake mackerels and cutlassfishes of the world (Families Gempylidae and Trichiuridae). Food and Agriculture Organization of the United Nations. 1993; 15.

2. Nelson JS, Grande TC, Wilson MVH. Fishes of the World. 5th ed. Hoboken, New Jersey: John Wiley & Sons Inc.; 2016.

3. FAO. Trichiurus lepturus (Linnaeus, 1758). Species Fact Sheet. 2021. https://www.fao.org/fishery/species/2468/en.

4. Kim S, Rho H. The fluctuation of fishing grounds of hairtail, Trichiurus lepturus in the Cheju strait by the distribution of feeding biomass. J Korean Soc Fish Ocean Technol. 2002;38: 1–10.

5. Zhang QY, Hong WS, Chen SX. Stock changes and resource protection of the large yellow croaker (Larimichthys crocea) and ribbon fish (Trichiurus japonicus) in coastal waters of China. J Appl Oceanogr. 2017;36(3): 438–445.

6. Lin HC, Tsai CJ, Wang HY. Variation in global distribution, population structures, and demographic history for four Trichiurus cutlassfishes. PeerJ 9:e12639. pmid:35003932

7. Wang HY, Dong CA, Lin HC. DNA barcoding of fisheries catch to reveal composition and distribution of cutlassfishes along the Taiwan coast. Fish Res. 2017;187: 103–109.

8. Wang HY, Heino MP. Adaptive and plastic variation in growth and maturation of the cutlassfish Trichiurus japonicus in the subtropical Pacific Ocean. Fish Bull. 2018;116(2): 171–182.

9. Shih NT, Hsu KC, Ni IH. Age, growth and reproduction of cutlassfishes Trichiurus spp. in the southern East China Sea. J Appl Ichthyol. 2011;27(6): 1307–1315.

10. Liu Y, Cheng J, Chen Y. A spatial analysis of trophic composition: a case study of hairtail (Trichiurus japonicus) in the East China Sea. Hydrobiologia. 2009;632(1): 79–90.

11. Weersing K, Toonen RJ. Population genetics, larval dispersal, and connectivity in marine systems. Mar Ecol Prog Ser. 2009;393: 1–12.

12. Lin X, Yang J, Guo J, Zhang Z, Yin Y, Song X, et al. An asymmetric upwind flow, Yellow Sea Warm Current: 1. New observations in the western Yellow Sea. J Geophys Res. 2011;116(C04026).

13. Oh KH, Lee S, Song KM, Lie HJ, Kim YT. The temporal and spatial variability of the Yellow Sea Cold Water Mass in the southeastern Yellow Sea, 2009–2011. Acta Oceanol Sin. 2013;32(9): 1–10.

14. Yang Y, Li K, Du J, Liu Y, Liu L, Wang H, et al. Revealing the subsurface Yellow Sea Cold Water Mass from satellite data associated with Typhoon Muifa. J Geophys Res. 2019;124(10): 7135–7152.

15. Tseng RS, Jan S, Zhu J. Circulation and Water Mass in the Taiwan Strait. Regional Oceanography of the South China Sea. 2020. pp. 433–470.

16. Ichikawa H, Beardsley RC. The current system in the Yellow and East China Seas. J Oceanogr. 2002;58(1): 77–92.

17. Pan HJ, Chen MT, Kong D, Lin X, Wong KT, Tsai HL, et al. Surface ocean hydrographic changes in the Western Pacific marginal seas since the Early Holocene. Front Earth Sci. 2020;8(200).

18. Dong YW, Wang HS, Han GD, Ke CH, Zhan X, Nakano T, et al. The impact of Yangtze River discharge, ocean currents and historical events on the biogeographic pattern of Cellana toreuma along the China coast. PLoS One. 2012;7(4): e36178. pmid:22563446

19. Wang J, Tsang LM, Dong Y-W. Causations of phylogeographic barrier of some rocky shore species along the Chinese coastline. BMC Evol Biol. 2015;15(1): 114. pmid:26071894

20. Tzeng CH, Chiu TY, Chen CS, Wang HY, Chiu TS. The current population structure of the demersal hairtail (Trichiurus japonicus) in the western North Pacific was shaped by the Taiwan Strait, as revealed by mitochondrial DNA. Taiwania. 2016;61(4): 279–287.

21. Carvalho GR, Hauser L. Molecular genetics and the stock concept in fisheries. Rev Fish Biol Fisheries. 1995;4: 55–79.

22. Ward RD. Genetics in fisheries management. Hydrobiologia. 2000;420(1): 191–201.

23. Utter FM. Biochemical genetics and fishery management: an historical perspective. J Fish Biol. 1991;39: 1–20.

24. Moritz C. Defining ‘evolutionarily significant units’ for conservation. Trends Ecol Evol. 1994;9(10): 373–375. pmid:21236896

25. He L, Zhang A, Weese D, Li S, Li J, Zhang J. Demographic response of cutlassfish (Trichiurus japonicus and T. nanhaiensis) to fluctuating palaeo-climate and regional oceanographic conditions in the China seas. Sci Rep. 2014;4(6380): 1–10. pmid:25223336

26. Xiao Y, Ren G, Song N, Li J, Gao T. Genetic diversity and population structure of largehead hairtail, Trichiurus japonicus, based on mtDNA control region. Mitochondrial DNA A: DNA Mapp Seq Anal. 2014;25(6): 473–481. pmid:23859051

27. Jarne P, Lagoda PJL. Microsatellites, from molecules to populations and back. Trends Ecol Evol. 1996;11(10): 424–429. pmid:21237902

28. Levinson G, Gutman GA. Slipped-strand mispairing: a major mechanism for DNA sequence evolution. Mol Biol Evol. 1987;4(3): 203–221. pmid:3328815

29. Tautz D. Hypervariability of simple sequences as a general source for polymorphic DNA markers. Nucleic Acid Res. 1989;17(16): 6463–6471. pmid:2780284

30. Sia EA, Jinks-Robertson S, Petes TD. Genetic control of microsatellite stability. Mutat Res. 1997;383(1): 61–70. pmid:9042420

31. Abdul-Muneer PM. Application of microsatellite markers in conservation genetics and fisheries management: recent advances in population structure analysis and conservation strategies. Genet Res Int. 2014;2014(691759). pmid:24808959

32. Zhan L, Paterson IG, Fraser BA, Watson B, Bradbury IR, Nadukkalam Ravindran P, et al. MEGASAT: automated inference of microsatellite genotypes from sequence data. Mol Ecol Resour. 2017;17(2): 247–256. pmid:27333119

33. Mardis ER. Next-generation DNA sequencing methods. Annu Rev Genomics Hum Genet. 2008;9(1): 387–402. pmid:18576944

34. Mikheyev AS, Vo T, Wee B, Singer MC, Parmesan C. Rapid microsatellite isolation from a butterfly by de novo transcriptome sequencing: performance and a comparison with AFLP-derived distances. PLoS One. 2010;5(6): e11212. pmid:20585453

35. Vukosavljev M, Esselink GD, van ‘t Westende WPC, Cox P, Visser RGF, Arens P, et al. Efficient development of highly polymorphic microsatellite markers based on polymorphic repeats in transcriptome sequences of multiple individuals. Mol Ecol Resour. 2015;15(1): 17–27. pmid:24893879

36. Postolache D, Leonarduzzi C, Piotti A, Spanu I, Roig A, Fady B, et al. Transcriptome versus genomic microsatellite markers: highly informative multiplexes for genotyping Abies alba Mill. and congeneric species. Plant Mol Biol Rep. 2014;32(3): 750–760.

37. Zhou R, Wu Z, Jiang FL, Liang M. Comparison of gSSR and EST-SSR markers for analyzing genetic variability among tomato cultivars (Solanum lycopersicum L.). Genet Mol Res. 2015;14(4): 13184–13194. pmid:26535631

38. Varshney RK, Graner A, Sorrells ME. Genic microsatellite markers in plants: features and applications. Trends Biotechnol. 2005;23(1): 48–55. pmid:15629858

39. Kim KS, Ratcliffe ST, French BW, Liu L, Sappington TW. Utility of EST-derived SSRs as population genetics markers in a beetle. J Hered. 2008;99(2): 112–124. pmid:18222933

40. González P, Dettleff P, Valenzuela C, Estrada JM, Valdés JA, Meneses C, et al. Evaluating the genetic structure of wild and commercial red cusk-eel (Genypterus chilensis) populations through the development of novel microsatellite markers from a reference transcriptome. Molecular biology reports. 2019;46(6): 5875–5882. pmid:31598817

41. Premachandra HKA, Lafarga-De la Cruz F, Takeuchi Y, Miller A, Fielder S, O’Connor W, et al. Genomic DNA variation confirmed Seriola lalandi comprises three different populations in the Pacific, but with recent divergence. Sci Rep. 2017;7(1): 9386. pmid:28839140

42. Wirgin I, Maceda L, Tozer M, Stabile JE, Waldman J. Atlantic coastwide population structure of striped bass Morone saxatilis using microsatellite DNA analysis. Fish Res. 2020;226: 105506. pmid:34987272

43. Cheng J, Zhang N, Sha Z. Nuclear microsatellites reveal population genetic structuring and fine-scale pattern of hybridization in the Japanese mantis shrimp Oratosquilla oratoria. PeerJ. 2020;8: e10270. pmid:33194430

44. Zhao Y, Zhu X, Li Z, Xu W, Dong J, Wei H, et al. Genetic diversity and structure of Chinese grass shrimp, Palaemonetes sinensis, inferred from transcriptome-derived microsatellite markers. BMC Genet. 2019;20(1): 75. pmid:31604423

45. Bi JZ, Shao CW, Miao GD, Ma HY, Chen SL. Isolation and characterization of 12 microsatellite loci from cutlassfish (Trichiurus haumela). Conserv Genet. 2009;10(4): 1171–1173.

46. An HS, Lee JH, Noh JK, Kim HC, Park CJ, Min BH, et al. Ten new microsatellite markers in cutlassfish Trichiurus lepturus derived from an enriched genomic library. Anim Cells Syst. 2010;14(3): 169–174.

47. Guo YS, Liu XM, Wang ZD, Lu HS, Liu CW. Isolation and characterization of microsatellite DNA loci from Naihai cutlassfish (Trichiurus nanhaiensis). J Genet. 2012;91(3): 109–112.

48. Liu S, Zu D, Liu Q, Dai F, Ma Q, Zhuang Z. Isolation and characterization of polymorphic microsatellite markers for Eupleurogrammus muticus. Conserv Genet Resour. 2015;7(2): 487–488.

49. Zhang H, Liang Z, Wu R, Niu S, Liang Y, Wang Q, et al. Microsatellite loci isolation in the Savalai hairtail (Lepturacanthus savala) based on SLAF-seq technology and generality in the related species. Genomics Appl Biol. 2018;37(8): 3331–3338.

50. Gu S, Yi MR, He XB, Lin PS, Liu WH, Luo ZS, et al. Genetic diversity and population structure of cutlassfish (Lepturacanthus savala) along the coast of mainland China, as inferred by mitochondrial and microsatellite DNA markers. Reg Stud Mar Sci. 2021;43: 101702.

51. Andrews S. FastQC—A Quality Control Tool for High Throughput Sequence Data. v0.11.5 ed. Babraham Institute, Cambridge, United Kingdom. 2010.

52. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15): 2114–2120. pmid:24695404

53. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7): 644. pmid:21572440

54. Haas BJ, Papanicolaou A, Yassour M, Grabherr MG, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-Seq: reference generation and analysis with Trinity. Nat Protoc. 2013;8(8): 1494–1512. pmid:23845962

55. Davidson NM, Oshlack A. Corset: enabling differential gene expression analysis for de novo assembled transcriptomes. Genome Biol. 2014;15(7): 410. pmid:25063469

56. Beier S, Thiel T, Münch T, Scholz U, Mascher M. MISA-web: a web server for microsatellite prediction. Bioinformatics. 2017;33(16): 2583–2585. pmid:28398459

57. Thiel T, Michalek W, Varshney R, Graner A. Exploiting EST databases for the development and characterization of gene-derived SSR-markers in barley (Hordeum vulgare L.). Theor Appl Genet. 2003;106: 411–422. pmid:12589540

58. Bachtrog D, Agis M, Imhof M, Schlötterer C. Microsatellite variability differs between dinucleotide repeat motifs—evidence from Drosophila melanogaster. Mol Biol Evol. 2000;17(9): 1277–1285. pmid:10958844

59. Dieringer D, Schlötterer C. Two distinct modes of microsatellite mutation processes: evidence from the complete genomic sequences of nine species. Genome Res. 2003;13(10): 2242–2251. pmid:14525926

60. Ashley MV, Dow BD. The use of microsatellite analysis in population biology: background, methods and potential applications. In: Schierwater B, Streit B, Wagner GP, DeSalle R, editors. Molecular ecology and evolution: approaches and applications. 69. Basel, Switzerland: Birkhäuser Verlag 1994. pp. 185–201.

61. Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, et al. Primer3—new capabilities and interfaces. Nucleic Acid Res. 2012;40(15): 115. pmid:22730293

62. Daniels J, Holmans P, Williams N, Turic D, McGuffin P, Plomin R, et al. A simple method for analyzing microsatellite allele image patterns generated from DNA pools and its application to allelic association studies. Am J Hum Genet. 1998;62(5): 1189–1197. pmid:9545387

63. Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004;4(3): 535–538.

64. Louis EJ, Dempster ER. An exact test for Hardy-Weinberg and multiple alleles. Biometrics. 1987;43(4): 805–811. pmid:3427165

65. Raymond M, Rousset F. GENEPOP (Version 1.2): Population genetics software for exact tests and ecumenicism. J Hered. 1995;86(3): 248–249.

66. Rousset F. Genepop’007: a complete re-implementation of the genepop software for Windows and Linux. Mol Ecol Resour. 2008;8(1): 103–106. pmid:21585727

67. Guo SW, Thompson EA. A Monte Carlo method for combined segregation and linkage analysis. Am J Hum Genet. 1992;51(5): 1111. pmid:1415253

68. Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10(3): 564–567. pmid:21565059

69. Beaumont MA, Nichols RA. Evaluating loci for use in the genetic analysis of population structure. Proc R Soc B: Biol Sci. 1996;263(1377): 1619–1626.

70. Excoffier L, Hofer T, Foll M. Detecting loci under selection in a hierarchically structured population. Heredity. 2009;103(4): 285–298. pmid:19623208

71. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984;38(6): 1358–1370. pmid:28563791

72. Belkhir K, Borsa P, Chikhi L, Raufaste N, Bonhomme F. GENETIX 4.05, logiciel sous Windows TM pour la génétique des populations. Université de Montpellier II, Montpellier, France. Laboratoire Génome, Populations, Interactions; 1996–2004.

73. Nei M. Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics. 1978;89(3): 583–590. pmid:17248844

74. Nei M, Tajima F, Tateno Y. Accuracy of estimated phylogenetic trees from molecular data. J Mol Evol. 1983;19(2): 153–170. pmid:6571220

75. Takezaki N, Nei M, Tamura K. POPTREEW: web version of POPTREE for constructing population trees from allele frequency data and computing some other quantities. Mol Biol Evol. 2014;31(6): 1622–1624. pmid:24603277

76. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2): 945–959. pmid:10835412

77. Earl DA, vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4(2): 359–361.

78. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14(8): 2611–2620. pmid:15969739

79. Peakall R, Smouse PE. GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes. 2006;6(1): 288–295.

80. Smouse PE, Peakall R. Spatial autocorrelation analysis of individual multiallele and multilocus genetic structure. Heredity. 1999;82(5): 561–573. pmid:10383677

81. Beerli P, Felsenstein J. Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. PNAS. 2001;98(8): 4563–4568. pmid:11287657

82. Dib C, Fauré S, Fizames C, Samson D, Drouot N, Vignal A, et al. A comprehensive genetic map of the human genome based on 5,264 microsatellites. Nature. 1996;380(6570): 152–154. pmid:8600387

83. Allendorf FW. Genetic drift and the loss of alleles versus heterozygosity. Zoo Biol. 1986;5(2): 181–190.

84. Lande R, Barrowclough GF. Effective population size, genetic variation, and their use in population. In: Soulé ME, editor. Viable populations for conservation. Cambridge University Press; 1987. pp. 87–123.

85. Cornuet JM, Luikart G. Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics. 1996;144(4): 2001–2014. pmid:8978083

86. Piry S, Luikart G, Cornuet JM. BOTTLENECK: a program for detecting recent effective population size reductions from allele data frequencies. J Hered. 1999;90: 502–503.

87. Luikart G, Allendorf FW, Cornuet JM, Sherwin WB. Distortion of allele frequency distributions provides a test for recent population bottlenecks. J Hered. 1998;89(3): 238–247. pmid:9656466

88. Kimura M, Crow JF. The number of alleles that can be maintained in a finite population. Genetics. 1964;49(4): 725. pmid:14156929

89. Di Rienzo A, Peterson AC, Garza JC, Valdes AM, Slatkin M, Freimer NB. Mutational processes of simple-sequence repeat loci in human populations. PNAS. 1994;91(8): 3166–3170. pmid:8159720

90. Shriver MD, Jin L, Chakraborty R, Boerwinkle E. VNTR allele frequency distributions under the stepwise mutation model: a computer simulation approach. Genetics. 1993;134(3): 983–993. pmid:8349120

91. Ockendon NF, O’Connell LA, Bush SJ, Monzón-Sandoval J, Barnes H, Székely T, et al. Optimization of next-generation sequencing transcriptome annotation for species lacking sequenced genomes. Mol Ecol Resour. 2016;16(2): 446–458. pmid:26358618

92. Hauser L, Seeb JE. Advances in molecular technology and their impact on fisheries genetics. Fish and Fisheries. 2008;9(4): 473–486.

93. Zhenzhen X, Ling X, Dengdong W, Chao F, Qiongyu L, Zihao L, et al. Transcriptome analysis of the Trachinotus ovatus: identification of reproduction, growth and immune-related genes and microsatellite markers. PLoS One. 2014;9(10): e109419. pmid:25303650

94. Du M, Li N, Niu B, Liu Y, You D, Jiang D, et al. De novo transcriptome analysis of Bagarius yarrelli (Siluriformes: Sisoridae) and the search for potential SSR markers using RNA-Seq. PLoS One. 2018;13(2): e0190343. pmid:29425202

95. Fang D, Zhou Y, Duan J, Zhang M, Xu D, Liu K, et al. Screening potential SSR markers of the anadromous fish Coilia nasus by de novo transcriptome analysis using Illumina sequencing. Genet Mol Res. 2015;14(4): 14181–14188. pmid:26600475

96. Duret L. Neutral theory: the null hypothesis of molecular evolution. Nat Educ. 2008;1(1): 803–806.

97. Kimura M. The neutral theory of molecular evolution: a review of recent evidence. Jpn J Genet. 1991;66(4): 367–386. pmid:1954033

98. Furlan E, Stoklosa J, Griffiths J, Gust N, Ellis R, Huggins RM, et al. Small population size and extremely low levels of genetic diversity in island populations of the platypus, Ornithorhynchus anatinus. Ecol Evol. 2012;2(4): 844–857. pmid:22837830

99. Hartl DL, Clark AG. Principles of population genetics. Sunderland, MA: Sinauer Associates; 1997.

100. Carlsson J. Effects of microsatellite null alleles on assignment testing. J Hered. 2008;99(6): 616–623. pmid:18535000

101. Wahlund S. Zusammensetzung von populationen und korrelationserscheinungen vom standpunkt der vererbungslehre aus betrachtet. Hereditas. 1928;11(1): 65–106.

102. Waples RS. Testing for Hardy–Weinberg proportions: have we lost the plot? J Hered. 2015;106(1): 1–19. pmid:25425676

103. Horikawa H, Zheng YJ, Meng TX. Biological and ecological characteristics of valuable fisheries resources from the East China Sea and the Yellow Sea. Seikai National Fisheries Research Institute, Japan. 2001.

104. De Meeûs T. Revisiting FIS, FST, Wahlund effects, and null alleles. J Hered. 2018;109(4): 446–456. pmid:29165594

105. Lewontin RC, Krakauer J. Distribution of gene frequency as a test of the theory of the selective neutrality of polymorphisms. Genetics. 1973;74(1): 175–195. pmid:4711903

106. Sun DQ, Li HY, Xu TJ, Wang RX. Development of microsatellite markers for the small yellow croaker Larimichthys polyactis (Sciaenidae) by cross-species amplification. Genet Mol Res. 2012;11(2): 1469–1474. pmid:22653596

107. DeWoody JA, Avise JC. Microsatellite variation in marine, freshwater and anadromous fishes compared with other animals. J Fish Biol. 2000;56(3): 461–473.

108. Xu D, Lou B, Zhou W, Chen R, Zhan W, Liu F. Genetic diversity and population differentiation in the yellow drum Nibea albiflora along the coast of the China Sea. Mar Biol Res,. 2017;13(4): 456–462.

109. Yang WT, Feng F, Yue GH. Isolation and characterization of microsatellites from a marine foodfish species ribbonfish Trichiurus haumela. Mol Ecol Notes. 2007;7(5): 781–783.

110. Kimura M. The neutral theory of molecular evolution. Sci Am. 1979;241(5): 98–129. pmid:504979

111. Avise JC. The history and purview of phylogeography: a personal reflection. Mol Ecol. 1998;7(4): 371–379.

112. Zhao X. In situ target-strength measurement of young hairtail (Trichiurus haumela) in the Yellow Sea. ICES J Mar Sci. 2006;63(1): 46–51.

113. Hsu KC, Shih NT, Ni IH, Shao KT. Genetic variation in Trichiurus lepturus (Perciformes: Trichiuridae) in waters off taiwan: several species or cohort contribution. Raff Bull Zool. 2007;14: 215–220.

114. Flowers JM, Schroeter SC, Burton RS. The recruitment sweepstakes has many winners: Genetic evidence from the sea urchin Strongylocentrotus purpuratus. Evolution. 2002;56(7): 1445–1453. pmid:12206244

115. Ott J, Rabinowitz D. The effect of marker heterozygosity on the power to detect linkage disequilibrium. Genetics. 1997;147(2): 927–930. pmid:9335624

116. Dias PC. Sources and sinks in population biology. Trends Ecol Evol. 1996;11(8): 326–330. pmid:21237863

117. Boughton DA. Empirical evidence for complex source-sink dynamics with alternative states in a butterfly metapopulation. Ecology. 1999;80(8): 2727–2739.

118. Guo C, Zhang G, Sun J, Leng X, Xu W, Wu C, et al. Seasonal responses of nutrient to hydrology and biology in the southern Yellow Sea. Cont Shelf Res. 2020;206: 104207.

119. Li Y, Wolanski E, Zhang H. What processes control the net currents through shallow straits? A review with application to the Bohai Strait, China. Estuar Coast Shelf Sci. 2015;158: 1–11.

120. Wang BD, Wang XL, Zhan R. Nutrient conditions in the Yellow Sea and the East China Sea. Estuar Coast Shelf Sci. 2003;58(1): 127–136.

121. Zheng LW, Zhai WD. Excess nitrogen in the Bohai and Yellow seas, China: Distribution, trends, and source apportionment. Sci Total Environ. 2021;794: 148702. pmid:34214818

122. Grant WAS, Bowen BW. Shallow population histories in deep evolutionary lineages of marine fishes: insights from sardines and anchovies and lessons for conservation. J Hered. 1998;89(5): 415–426.

123. Hewitt G. The genetic legacy of the Quaternary ice ages. Nature. 2000;405(6789): 907–913. pmid:10879524

124. Lambeck K, Esat TM, Potter EK. Links between climate and sea levels for the past three million years. Nature. 2002;419(6903): 199–206. pmid:12226674

125. Avise JC. Phylogeography: the history and formation of species: Harvard university press; 2000.

126. Liu JX, Gao TX, Wu SF, Zhang YP. Pleistocene isolation in the Northwestern Pacific marginal seas and limited dispersal in a marine fish, Chelon haematocheilus (Temminck & Schlegel, 1845). Mol Ecol. 2007;16(2): 275–288. pmid:17217344

127. Schulte PM. Environmental adaptations as windows on molecular evolution. Comp Biochem Physiol B, Biochem Mol Biol. 2001;128(3): 597–611. pmid:11250554

128. Han Z, Yanagimoto T, Zhang Y, Gao T. Phylogeography study of Ammodytes personatus in Northwestern Pacific: Pleistocene isolation, temperature and current conducted secondary contact. PLoS One. 2012;7(5): e37425. pmid:22685527

129. Wang P. Response of Western Pacific marginal seas to glacial cycles: paleoceanographic and sedimentological features. Mar Geol. 1999;156(1–4): 5–39.

130. Ni G, Li Q, Kong L, Yu H. Comparative phylogeography in marginal seas of the northwestern Pacific. Mol Ecol. 2014;23(3): 534–548. pmid:24600706

131. Zhang K, Li J, Hou G, Huang Z, Shi D, Chen Z, et al. Length-based assessment of fish stocks in a data-poor, jointly exploited (China and Vietnam) fishing ground, northern South China Sea. Front Mar Sci. 2021;8: 1043.

132. Ding Q, Shan X, Jin X, Gorfine H. A multidimensional analysis of marine capture fisheries in China’s coastal provinces. Fish Sci. 2021;87(3): 297–309.

133. Makino M, Watari S, Hirose T, Oda K, Hirota M, Takei A, et al. A transdisciplinary research of coastal fisheries co-management: the case of the hairtail Trichiurus japonicus trolling line fishery around the Bungo Channel, Japan. Fish Sci. 2017;83(6): 853–864.

134. Watari S, Tokumitsu S, Hirose T, Ogawa M, Makino M. Stock structure and resource management of hairtail Trichiurus japonicus based on seasonal broods around the Bungo Channel, Japan. Fish Sci. 2017;83(6): 865–878.

135. Ling JZ, Li SF, Yan LP, Cheng JH. Utilization and management of Trichiurus japonicus resources in East China Sea based on Beverton-Holt model. Ying Yong Sheng Tai Xue Bao. 2008;19(1): 178–182. pmid:18419092

136. Padhi BK, Mandal RK. Applied fish genetics. India, Visakhapatnam: Fishing Chimes; 2000.