A total of 135 P. echinogaster samples (Figure 1) were collected at nine different locations along the coast of China (15 individuals at each location) from October 2017 to December 2017. External morphological identification of the samples was performed mainly by referring to Nakabo (2013) and Li et al. (2019). For each fresh P. echinogaster sample, sterile scissors and tweezers were used to obtain back muscle tissues. All muscle tissues were stored in a cryogenic freezer at -80°C for future experiments. The phenol–chloroform method was used to extract genomic DNA from each sample (Li et al., 2019; Xue et al., 2019). Then, 1% agarose gel electrophoresis and an Invitrogen Qubit fluorometer were used to assess the degradation degree and concentration of the genomic DNA. Qualified genomic DNA (500 ng) was submitted to Biomarker Technologies Corporation (Beijing, China) for library construction and sequencing.

According to the genome size and guanine+cytosine (G+C) content of P. argenteus, enzymatic digestion prediction was performed. The analysis software for SLAF enzymatic digestion prediction independently developed by Biomarker Technologies Corporation (Sun et al., 2013) was used to predict the digestion of the reference genome, and the optimal digestion scheme was selected according to the following principles: (1) the percentage of enzymatic fragments located in the repeat sequences was as low as possible; (2) the enzymatic fragments were distributed as evenly as possible across the genome; (3) the length of each enzymatic fragment was consistent with that of the specific experimental system; and (4) the number of obtained enzymatic fragments (SLAF tags) agreed with the expected number of tags. With the HaeIII restriction enzyme, the genomic DNA of each qualified sample was digested separately. Each obtained digestion fragment (SLAF tag) had a poly (A) tag added to the 3’ end, and dual-index sequencing adaptors were attached (Kozich et al., 2013). Genomic DNA was incubated at 37°C with MseI, T4 DNA ligase, ATP, and MseI adapter. These restriction-ligation reactions for further PCRs were heat-inactivated at 65°C and then digested in an additional reaction with the restriction enzymes HaeIII and BfaI at 37°C. Furthermore, PCR was performed with diluted restriction-ligation PCR products, dNTP, Taq DNA polymerase and MseI primers containing barcode1. The further PCR products were purified using an E.Z.N.A Cycle Pure Kit. The purified PCR products were pooled and incubated at 37°C with MseI, T4 DNA ligase, ATP, and Solexa adapter. The reaction products after incubation were purified using a Quick Spin column and then electrophoresed on a 2% agarose gel. A gel extraction kit was applied to select SLAFs 314-394 bp in size. The gel-purified fragments were further amplified by PCR with Phusion Master Mix and Solexa Amplification primer mix to add barcode2. Finally, paired-end sequencing was performed on the Illumina HiSeq 2500 platform. To evaluate the accuracy of enzyme digestion, Nipponbare was selected as a control for sequencing.

The raw data obtained from sequencing were characterized using the dual-index sequencing adaptor. The raw reads of each sample were obtained, and reads with a quality score lower than 30 were excluded. Based on sequence similarity, the remaining high-quality reads of each sample were clustered, and the reads with a similarity greater than 98% were considered to have clustered as a single SLAF tag. A SLAF tag with sequence differences between samples was defined as polymorphic. The sequence with the greatest sequencing depth for each SLAF tag was taken as the reference sequence. The Burrows–Wheeler alignment tool was used (Li and Durbin, 2009) to align the reads to the SLAF tags, both GATK (parameters: –stand_call_conf 30.0 and –stand_emit_conf 0) (McKenna et al., 2010) and SAMTOOLS (parameters: default) (Li et al., 2009) were used to perform single nucleotide polymorphism (SNP) calling, and the overlapping SNPs in the GATK and SAMTOOLS results were used as the final SNP dataset. The generated SNPs were saved in a variant call format (VCF) file. To ensure the accuracy of subsequent analyses, VCFtools (Danecek et al., 2011) was used to screen SNPs with the following parameters: -MAF 0.01 (minimum allele frequency> 0.01); –max-missing 0.1 (filtering out the genotypes with less than 90% data); –min-meanDP 150 (minimum mean depth of coverage > 150); –min-alleles 2 –max-alleles 2 (only two alleles); –minGQ 98 (quality score> 98); –minQ 30 (retaining the loci with a quality score> 30); –remove-indels (excluding the loci containing indels); and -HWE 0.05 (excluding the loci with P < 0.05 in the Hardy-Weinberg equilibrium test).

To describe the genetic diversity levels of all genome-wide SNPs of the nine P. echinogaster locations, TASSEL software (version 5.2.31) (Bradbury et al., 2007) was used to estimate the nucleotide diversity (Pi) and Tajima’s D value of each SNP in each P. echinogaster population, the step size and window size was 10 and 20, respectively. With Circos software (Krzywinski et al., 2009), the Pi, expected heterozygosity (H_E), and Tajima’s D of each SNP were visualized.

Using the “populations” module of Stacks software (version 1.34) (Catchen et al., 2013), the genetic diversity levels of the nine populations, including the polymorphic loci Pi, observed heterozygosity (H_O), H_E, and inbreeding coefficient (F_IS), were statistically analyzed. Arlequin software (version 3.5.2.2) (Excoffier and Lischer, 2010) was used to estimate the pairwise genetic differentiation (F_ST) of P. echinogaster, and 10,000 permutations were used to test the significance.Based on all SNPs, four methods were used to estimate population structure and individual clustering within a population. (1) PGDSpider software (version 2.0.5.2) (Lischer and Excoffier, 2012) was used to convert the VCF file of all SNPs to STR format, and then ADMIXTURE software, based on the maximum likelihood method (version 1.3.0) (Alexander et al., 2009), was used to evaluate ancestors. This software uses a fast numerical optimization algorithm to achieve a fast analysis rate. During analysis, the range of genetic clusters (K) was set to 1-9, and each analysis was repeated 10 times. Finally, based on the cross-validation (CV) error corresponding to the optimal K value, the clustering results were plotted. (2) Using TASSEL software (version 5.2.31) (Bradbury et al., 2007), neighbour-joining (NJ) trees of different P. echinogaster geographical populations were constructed to clarify the phylogenetic relationships between all individuals in the P. echinogaster geographical populations, and the number of different residues between sequences was used to measure genetic distance. The acquired NJ trees were visualized with iTOL software (https://itol.embl.de/). (3) The adegenet package of R software (Jombart et al., 2010) was used to perform principal component analysis (PCA) of all individuals and visualize the geographical populations and interpopulation relationships. (4) NetView P software, which adopts the k-nearest neighbour (KNN) algorithm, was used to visualize the network topology among all individuals to accurately and thoroughly explore the refined population structure of P. echinogaster. Before running NetView P software, the range of KNN values was first set to 1-60, and the optimal KNN value was determined according to the Fast-Greedy, Infomap, and Walktrap algorithms. In other words, the genetic similarity of individuals was determined based on the optimal resolution of the genetic structure, and then the network topology of all individuals was obtained based on the optimal KNN value. (5) Arlequin (version 3.5.2.2) (Excoffier and Lischer, 2010) was used for analysis of molecular variance (AMOVA) to detect differentiation between groups (F_CT) and differentiation between geographical populations within a group (F_SC), and the number of iterations was set to 5000. In this study, the P. echinogaster geographical populations were divided into three groups based on geographic location: Dalian, Tianjin, Qingdao, and Nantong; Zhoushan and Xiamen; and Shantou, Zhuhai, and Zhanjiang. Additionally, the effective population size (Ne) of nine P. echinogaster geographical populations was estimated with NeEstimator v2 software (version 2) (Do et al.,  2014).

PGDSpider (version 2.0.5.2) (Lischer and Excoffier, 2012) was used to convert the VCF file of all SNPs to the GENEPOP format, and the gene flow between the nine geographical populations was analyzed using divMigrate-Online (https://popgen.shinyapps.io/divMigrate-online/). The parameters were as follows: the number of bootstraps was 1000, alpha was 0.05 and the migration statistic was Nm.

The genotype-environment association method implemented in Bayenv (version 2.0) (Coop et al., 2010) was applied to detect putative SNPs correlated with temperature variations. First, we obtained the high-resolution mean lowest sea surface temperature (LSST), mean annual sea surface temperature (ASST) and mean highest sea surface temperature (HSST) data of nine sea areas over 68 years (from 1950 to 2017) by combining data from the Japan Meteorological Agency (JMA; http://www.jma.go.jp/jma/index.html), Advanced Very High Resolution Radiometer (AVHRR; http://oceanwatch.pifsc.noaa.gov/thredds/catalog.html) and Geostationary Operational Environmental Satellites (GOES; http://oceanwatch.pifsc.noaa.gov/thredds/catalog.html). We first converted SNP files and average environmental data files (Table 1) into matrix files needed for software operation, and the matrix was evaluated according to the software specification. The parameters were set as follows: -k 100000 -r 63479. Then, Bayenv tests were applied to identify putative SNPs correlated with temperature variations using the “calc_bfs.sh” script in Bayenv2 software, and a Bayes factor (BF) value higher than 20 (Kass and Raftery, 1995) was set as the filtering condition for putative SNPs. Thereafter, we used the overlapping ASST-related, LSST-related, and HSST-related SNPs as candidate temperature-selected SNPs. SLAF sequences containing these SNPs were then annotated using Blast2GO software (Conesa et al., 2005).

After filtering, a total of 448.17 million high-quality paired-end reads were obtained from SLAF-seq of 135 P. echinogaster samples. The average Q30 of the reads was 95.59%, and the average GC content was 43.17%. In addition, for the control sequence of Nipponbare, which was used to evaluate the accuracy of the established library, 1.79 M reads were obtained. After clustering high-quality reads, we obtained a total of 935,215 SLAF tags, and the average sequencing depth of the SLAF tags was 20.80×. Among the SLAF tags, 461,932 were polymorphic. With GATK and SAMTOOLS, 6,152,196 SNPs were further obtained. After the removal of low-quality SNPs, we obtained 62,473 SNPs, and 10,707 passed the Hardy-Weinberg equilibrium test (P < 0.05). After retaining one SNP from each locus with the lowest missing rate, 5,579 SNPs were excluded. Finally, 46,187 highly consistent SNPs were identified for subsequent genetic structure analysis.

The Pi and Tajima’s D values of the 46,187 SNPs exhibited similar fluctuating patterns among the nine P. echinogaster geographical populations (Figure 2). The Tajima’s D values of most SNPs in the nine P. echinogaster geographical populations were less than 0, indicating that many rare alleles existed at a high frequency, which may have been due to the large population size.

After filtering, a total of 46,103 SNPs were considered neutral SNPs, meaning those SNPs that selected balance and genotype were approximately constant. Statistical results of genetic diversity at the 46,103 SNPs showed that the average H_O, average H_E, and Pi varied among the nine P. echinogaster locations (H_O = 0.2046-0.2364, H_E = 0.2222-0.2337, and Pi = 0.2336-0.2425), and the F_IS values of all nine geographical populations were small (F_IS = 0.0222-0.0806) (Table 2). The Shantou population had the lowest average genetic diversity (H_O, H_E, and Pi) and Tianjin showed the highest genetic diversity (Table 2). The percentage of polymorphic SNPs in each of the nine P. echinogaster geographical populations was relatively high (86.5179-93.2759%).

The nine P. echinogaster geographical populations had low values of genetic differentiation, ranging from -0.00109 to 0.00267 (Table 3). ADMIXTURE analysis of 46,103 SNPs revealed an optimum K value of 1 with a minimum CV error. The clustering results from ADMIXTURE based on K=1 showed that all the individuals in the nine geographical locations were clustered into one group, and the clustering results from ADMIXTURE based on K = 2 to 9 showed a similar result (Figure 3A). The clustering pattern from ADMIXTURE was validated by analysis of NJ trees (Figure 3B) and the PCA results (Figure 3C) for all 46,103 SNPs. In addition, NetView P software revealed fine genetic structure among the nine P. echinogaster geographical populations. In this study, the optimal KNN value was determined to be 38 based on various algorithms, and the network topological results showed that all the individuals in all nine locations were clustered together (Figure 3D). The AMOVA results (Table 4) also indicated that most of the variation was within populations (F_ST = 0.03788; P = 0.00098), while statistically nonsignificant genetic differentiation was observed among the three groups (Va = 0.01; F_CT = 0.00014, P = 0.12317; Table 4).

divMigrate-Online software was used to analyse the gene exchange between the nine P. echinogaster geographical populations, and the results showed extensive gene exchange (0.34 - 1) among them (Figure 4), providing good support for the above results on the lack of genetic differentiation.

In the present study, we first calculated the average sea surface temperature (ASST), low sea surface temperature (LSST) and high sea surface temperature (HSST) of nine sea areas over 68 years (Table 1). Then, Bayenv software revealed a total of 329 SNPs associated with temperature variables. Of these SNPs, 203 were associated with ASST, 190 were associated with LSST, and 202 were associated with HSST. There was 84 overlapping SNPs were associated with ASST, LSST, and HSST and we used the 84 overlapping ASST-related, LSST-related, and HSST-related SNPs as candidate temperature-selected SNPs (Figure 5; Table 3). The genetic differentiation of nine P. echinogaster geographical populations based on temperature-selected SNPs was also low (Figure 6). Whole-genome sequences containing 84 SNPs were then used for further annotations, and the results showed that 51 gene sequences containing temperature-selected SNPs matched homologous protein sequences in the nonredundant protein sequences (Nr) database, and these genes were primarily related to material metabolism, ion transport, breeding, stress reactions, and inflammatory reactions. Next, the enrichment of sequences containing temperature-selected SNPs in Gene Ontology (GO) categories was tested, and 30 significantly enriched GO terms were identified (Figure 7).

This is the first study exploring the genetic structure and the genetic diversity of different geographical populations of P. echinogaster at the genomic level. Although the geographic populations of P. echinogaster had higher H_O, H_E, and Pi values than those of Pampus minor, which we previously studied (Li et al., 2021), they were still lower than 0.50. Except for the ST population, the F_IS values of the other 8 geographical populations were all lower than the standard maximum value (0.05) of no genetic differentiation between geographical populations. This means that the genetic differentiation degree between the ST population and the other 8 geographical populations is relatively large.

Although P. echinogaster occupies a wide variety of habitats, this study confirmed that there was no significant genetic differentiation between P. echinogaster geographical populations based on multiple analytical methods. Gene flow analysis also showed frequent gene exchange between P. echinogaster geographical populations, which can be explained by random mating and suggests that high genetic homogeneity exists between P. echinogaster geographical populations. These results are consistent with our previous results based on mitochondrial DNA sequences and microsatellite loci (Li et al., 2018). We speculate that the population expansion of this species after the last glacial period and the lack of a balance between migration and genetic drift may be two main reasons for their unclear genetic structure (Wu et al., 2012; Li et al., 2018). Similar genetic structures are also found in species with closed expansion times, such as Trachurus japonicus (Song et al., 2013), Sebastes schlegelii (Zhang et al., 2013), Synechogobius ommaturus (Song et al., 2010), and Larimichthys polyactis (Xiao et al., 2009; Xiao et al., 2013). At the same time, the open marine environment lacks obvious physical barriers; the strong diffusion ability of eggs, larvae, and adults; and the high randomness of diffusivity provided by ocean currents are key reasons for the genetic homogeneity of marine organisms (Hewitt, 2000). P. echinogaster is a seasonal migratory fish species that lays eggs in batches. The spawning period can last for more than two months, and the number of eggs per mother can be 117,000-218,000 (Jin et al., 2005). In addition, the eggs of P. echinogaster are buoyant, and the larvae have a pelagic period of at least 1 month. These life history characteristics may be conducive to the long-distance diffusion of eggs and larvae and can promote gene exchange between P. echinogaster geographical populations, ultimately resulting in very low population genetic differentiation within a very wide distribution range (Grant and Bowen, 1998). These ocean currents can also affect the dispersal of marine larvae. For instance, the cold and low-salinity China Coastal Current (CCC) flows southwards along the coast of China from the Bohai Sea to the Taiwan Strait (Wu, 1982). The South China Sea Warm Current (SCSWC) connects the East China Sea and the South China Sea (Wu, 1982). Additionally, the Kuroshio Current transports warm and high salinity water northeastward from the Lozon Strait to the steep east-west continental shelf rupture in the East China Sea (Cheng and Sha, 2017). In conclusion, we speculate that the complex ocean current system along the coast of China may further affect the population size and genetic connectivity of this species by increasing the diffusion of eggs and larvae. Currently, P. echinogaster geographical populations form a continuous distribution and present a population-wide genetic pattern (Li et al., 2018). In summary, the life history characteristics of P. echinogaster and the influence of marine environmental factors may cause the mixing of P. echinogaster geographical populations, thus resulting in frequent gene exchange and high genetic homogeneity between them.

Temperature-selected SNPs can be applied to explore the prospects of P. echinogaster individual allocation to their region of origin (Nielsen et al., 2009; Jasonowicz et al., 2016). Considering that P. echinogaster is a seasonal migratory fish with long-distance diffusion (Yamada et al., 2007; Li et al., 2018), aquatic temperature conditions vary among geographical populations of this species. In fact, we found that there were differences in the high-resolution mean annual sea surface temperature (ASST), highest sea surface temperature (HSST), and lowest sea surface temperature (LSST) of nine sampling sites in 68 years (1950–2017). The results showed that multiple temperature-selected SNPs were related to material metabolism. In fact, the metabolic capacity of fishes is particularly susceptible to environmental stress; thus, fish use critical energy materials (amino acids, fat, and sugar) stored in their bodies to adapt to changes in habitat factors and for the energy required for life activities, such as migration (Pastoureaud, 1991; Wang, 1999; Liang, 2005). On the other hand, fish also regulate the fatty acid composition of their bodies to increase membrane fluidity, thereby adapting to changes in habitat temperature (Li et al., 2019). Therefore, genes related to material metabolism may play an important role in the adaptation of P. echinogaster to the temperatures of different geographical environments. Ion (such as calcium) transport signalling pathways also play important roles in the process of local adaptive evolution in P. echinogaster. For example, calcium ions are involved in fish reproduction, development, learning and memory, mitochondrial function, muscle contraction, and other functions (Berridge et al., 2000). Calcium signalling is also thought to play an important role in the regulation of ion exchange and osmoregulation, which makes this pathway a reasonable target of spatially differentiated selection under osmotic pressure in P. echinogaster because this species feeds on a wide variety of prey. In fact, a previous study confirmed that temperature changes the osmotic pressure regulation of marine organisms by altering the salinity of seawater (Whiteley et al., 2018). However, these results may have been due to chance because we detected only a small number (n = 84) of temperature-selected SNPs. Therefore, we still need to capture as many highly divergent regions of the genome and temperature-selected SNPs as possible (Woodings et al., 2018). Moreover, we have to acknowledge that explicitly linking genomic features of local adaptation to specific environmental conditions is not an easy task, and further targeted research is needed (Woodings et al., 2018).