Results
We sampled 49 individuals from five populations covering the East China Sea, the Yellow Sea, and the Pacific coastal waters of Japan (Figure 1a and Table 1). Whole-genome resequencing yielded 680 Gb of sequencing data and generated 13.6 Gb for each individual. Genome alignment resulted in an average depth of 20.91-fold and an average of 99.40% genome coverage (at least one fold) (Supplementary Table S1). All individual data were aligned to the reference genome (mapping rate from 95.00% to 96.89%) (Supplementary Table S2), and SNPs were called after rigorous quality filtering. We identified 5,483,086 SNPs after quality control for further analysis, with 406,925 in exonic regions, 2,025,230 in intronic regions, and 2,649,227 in intergenic regions. Of the exonic SNPs, we identified 312,482 synonymous and 94,443 nonsynonymous SNPs (Supplementary Table S3). The nucleotide diversity π for each population was similar, ranging from 1.64 × 10–3 to 1.69 × 10–3 (Figure 1b).
Clustering analyses using PLINK was performed to examine the relationships among the populations (Figure 1c). The admixture analyses revealed clear evidence for clustering. With K = 2, the admixture structure showed two clusters. Populations from China (RS, ZS, and ST) formed one cluster, whereas populations in Japan (IB and TB) formed another cluster. When K = 3, the RS population formed a separate cluster. As the K value increased, the ST population separated as a unique cluster. When K increased from 3 to 6, the ZS population was shared with clusters from RS and ST. TB and IB were highly admixed in all cases with K from 2 to 6. PCA results recovered the similar clusters achieved by the admixture analyses (Figure 1d). The first principal component of PCA (PC1) separated China and Japan clades, consistent with the admixture result at K = 2. The second principal component of PCA (PC2) further separated RS, ZS, and ST populations but was not able to distinguish between TB and IB populations, consistent with the admixture results at K = 5 and 6 (Figure 1c). All individuals from China groups clustered within the population defined by their sampling locations, revealing a signal of geographic structure from the East China Sea to the Yellow Sea.
The phylogenetic tree further supported the admixture analysis and PCA results. A neighbor-joining tree was constructed on the basis of whole-genome SNPs (Figure 1e). The tree formed two clades that represent the geographic divergence between the coastal waters of China and Japan. In the China clade, three distinct clusters generally defined by geographic localities were recovered. However, the phylogenetic topology of the Japan clade showed a shallow structure, and some of the individuals were grouped incorrectly into another geographic population.
We calculated pairwise F ST between the populations to quantify their genetic differentiation (Table 2). Pairwise F ST ranged from −0.0001 to 0.0377, with an average of 0.0224, consistent with the overall structured population.F ST between China populations (RS, ZS, and ST) and Japan populations (IB and TB) ranged from 0.0237 to 0.0377, with an average of 0.0316. The level of genetic differentiation among the China populations was higher than that between the two Japan populations.
To explore the influence of geographic distance on genetic differentiation, we performed Mantel tests to test the association between the geographic distance matrix and the pairwiseF ST matrix. Two patterns of geographic distance (i.e., coastal distance and oceanic distance) were used in the Mantel tests (Figure 2). We identified a closer relationship (r = 0.90,P = 0.0002) betweenFST /(1−FST ) and coastal distance than that between oceanic distance andFST /(1−FST ) (r = 0.80,P = 0.0029), indicating isolation due to costal distance, with coastal distance explaining 91% of the variation in genetic differentiation for the species. These results of population structure analysis and Mantel tests collectively indicated that strong barriers might have had a greater influence on population differentiation than the other factors.
The demographic history of S. japonica was inferred to understand its evolutionary history. The historical effective population sizes of the populations were estimated via PSMC. Results showed that the demographic history of S. japonica could be traced back to approximately 100,000 years ago, and the China and Japan populations experienced evidently different demographic histories (Figure 3). Three China populations experienced similar demographic trajectories: one large population expansion and one population bottleneck. Three China populations suffered from an obvious population expansion approximately 30,000 years (Ka) ago. The effective population size of the China populations peaked ~2–10 Ka, and then the effective population sizes, especially of the RS population, declined ~2 Ka and formed a bottleneck. Among the three China populations, the signal of population expansion increased from north to south, but the signal of bottleneck decreased from north to south. The IB and TB Japan populations underwent a different demographic history compared with the China populations: only one small population expansion. In contrast to the large and rapid population expansion of the China populations, the Japan populations started to increase slowly ~10 Ka and reached the climax in the present.
We then reconstructed a maximum likelihood (ML) tree of the five populations by using TreeMix to address population history relationships and identify pairs of populations that are related to each other independent of that captured by this tree (Figure 4a). The ML without migration events inferred from the TreeMix analysis divided the 49 individuals into two clusters that were similar to the population structuring patterns identified from the phylogenetic tree analysis, PCA, and genetic structure analysis. When potential migration edges (i.e., migration events between the branches) were added to the ML tree, strong migration events were detected between the clusters (Figure 4b). We observed a statistically significant migration edge (P< 2.2E–308) with the estimated weight of 28.8%; this edge provided evidence for the gene flow from the Japan TB population into the RS population.
We compared the genomes of S. japonica individuals to identify signatures of positive selection under different environment selection pressures. Considering that genetic isolation occurred between the China and Japan populations, we chose the ST and ZS populations as the control groups to reveal candidate cold-temperature selection genes in the RS population. Using the top 5% maximum F ST(F ST ≥ 0.1260) and π ratio (ΠST/RS ≥ 1.2380) methods, a total of 132 candidate genes corresponding to 3.64 Mb in size were identified in the RS population, with the ST population as the control group (Figure 5a). Using the ZS population as the control group, 509 candidate genes were selected corresponding to 18.96 Mb in size (F ST ≥ 0.0904, ΠZS/RS ≥ 1.2676). To validate further the candidate genes under strong selective sweeps in the RS population, we recognized a total of 81 genes in the overlap of RS/ZS and RS/ST pairs as potentially affected genes associated with cold-temperature adaptation for subsequent analyses (Figure 5b, Supplementary Table S4).
With the low-temperature populations RS and ZS as the control groups and using the top 5% maximum F ST values (F ST ≥ 0.0485, 0.0845) and π ratio values (ΠRS/ST ≥ 1.2404, ΠZS/ST ≥ 1.1544), 1037 and 534 candidate genes associated under selective sweeps in ST population were identified, respectively. A total of 365 genes shared by the ST/RS and ST/ZS pairs associated with high-temperature adaptation in the ST population were identified (Figure 5c, Supplementary Table S5). Only two cold-temperature adaptation genes overlapped with the high-temperature adaptation genes. Overall, 81 and 365 genes associated with cold- and high-temperature environments were positively selected, respectively.
To detect possible parallel adaptation between ZS and IB/TB populations, we identified 692 genes in ZS and 122 genes in Japan groups involved in warm-temperature adaptation, with the cold-temperature population RS as the control group. A total of 111 out of 122 (91.0%) warm-temperature adaptation genes in the Japan populations overlapped with the warm-temperature adaptation genes of the ZS population (Supplementary Table S6), and no gene overlapped with the high-temperature adaptation genes of the ST population (Figure 5d). Considering the geographic isolation at the whole-genome SNPs and similar temperature environments, the highly shared selection genes between the Japan and ZS populations may be a possible evidence for parallel adaptive evolution.
The PCA results obtained on the basis of SNPs of the candidate genes related to cold-temperature adaptation showed that individuals from the RS population were correctly separated from the other populations (Figure 6a). PC1 tended to separate populations with different latitudes from south to north. PC2 separated the RS population from the other populations. Compared to high-temperature population, three warm-temperature populations (namely, ZS, IB, and TB) showed a closer relationship. Similarly, the PCA results that included the SNPs of the candidate genes for high-temperature adaptation separated the ST population (high temperature) from the other populations, with some individuals from the RS population at an intermediate position (Figure 6b). The allele frequency of one SNP within the cold-temperature adaptation gene Picalm  was significantly associated with population temperature (Figure 5e). However, the allele frequencies of one SNP within the warm-temperature adaptation gene SORCS3 showed high values in three warm-temperature populations (Figure 5f).
The PCA results that included the SNPs of the candidate genes from warm-temperature populations also divided the 49 individuals into two clades between the ST population and the other four populations, with the RS population at the intermediate position (Figure 6c). This topological pattern was different from that inferred from the genome-wide SNPs (Figure 1d), thus supporting the credibility of the candidate genes under selection. The substructure of the PCA based on the candidate parallel adaptation genes revealed a closer relationship between the ZS populations and IB/TB pair than with genome-wide SNPs, where the populations within the pair had large geographical and similar temperature environments (Figure 3, Table 1).
To obtain a broad overview of the molecular functions of these genes and detect the most differentiated regions of S. japonica genome, we performed GO enrichment and KEGG analyses on the candidate selection genes. These analyses offered a clear insight into the genetic evolution and adaptive mechanisms of S. japonica under different thermal environments.
The functional enrichment found that the candidate genes under cold selection were significantly enriched in only one KEGG pathway, that is, choline metabolism in cancer (ko05231, FDR adjusted P = 0.0421) (Figure 7, Supplementary Table S7). The GO term enrichment analysis for the 81 candidate genes under cold temperature selection showed that the genes were classified into 47 categories, including metabolic processes, biological regulation, response to the stimulus, and signaling process in the biological process; membrane and membrane part in the cellular component; and binding and catalytic activity in the molecular function (Supplementary Figure S1). We found 50 significantly enriched GO terms after FDR correction (Supplementary Figure S2, Supplementary Table S8). The significantly enriched GO categories were primarily associated with cell part, cation transmembrane transporter activity, tissue homeostasis, transport, and lipid metabolism.
In the ST population, 365 genes were identified under selection, and these genes were significantly enriched in the cell adhesion molecule pathway (ko04514, FDR adjusted P = 0.0046), tight junction pathway (ko04530, FDR adjusted P = 0.0189), and leukocyte transendothelial migration (ko04670, FDR adjusted P = 0.0189) (Figure 7, Supplementary Table S9). The selected genes were classified into 53 GO categories (Supplementary Figure S3), covering the major categories found in the RS population, including seven other categories (growth and behavior in the biological process; membrane-enclosed lumen in the cellular component; and transcription factor activity and protein binding in the molecular function). Four significantly enriched GO terms, namely, calcium-independent cell–cell adhesion via plasma membrane cell-adhesion molecules (GO: 0016338), establishment of the skin barrier (GO: 0061436), positive regulation of wound healing (GO: 0090303), and regulation of water loss via skin (GO: 0033561) were detected, all of which belong to the biological process (Supplementary Figure S2, Supplementary Table S10). Despite no significant GO terms in the cellular component and molecular function, the enriched categories in the cellular component and molecular function were mainly related to the cytoskeleton, the component of membrane, and transmembrane transporter activity. In this population, 13 genes in the cell-adhesion molecule pathway showed signals of natural selection. The 13 genes were nine members of Claudin family, JAM1 , CDH4 , andNLGN . Most members of the Claudin family were also located in the tight junction pathway and leukocyte transendothelial migration way. The enrichment of selected genes in the ST population suggested the importance of adaptation to high-temperature climate. The different patterns of enrichment for GO and KEGG pathways between the RS and ST populations may reflect the difference between cold- and warm-temperature stress.
A total of 111 genes were identified for parallel adaptations in the ZS, IB, and TB populations. These genes were classified into 48 categories (Supplementary Figure S4). Nineteen significantly overrepresented GO categories were detected for thermal adaptation (Supplementary Table S11). The GO clusters were primarily enriched in the categories of cell projection, structure of the cytoskeleton, binding, and maintenance of membrane (Supplementary Figure S2). Although no significantly enriched KEGG pathway was detected, the top 20 enriched pathways were related to metabolism (synthesis and degradation of ketone bodies, cyanoamino acid metabolism, linoleic acid metabolism, alpha-linolenic acid metabolism, and butanoate metabolism), circulatory system (vascular smooth muscle contraction), and endocrine system (salivary secretion and pancreatic secretion) (Figure 7, Supplementary Table S12). The enrichment of the selected genes in the warm-temperature population offered a parallel adaptive mechanism under thermal selection.
The receiver operating characteristic curve (AUC) results (0.983 in China group, 0.995 in Japan group) suggested that the Maxent model provides good predictive performance. The Maxent predictions for two groups of Japanese whiting suggested that the potential distribution of this species may greatly change (Figures 8 and 9). The predictions showed that a large part of coastal areas of China, central Japan, and Korea were suitable for the China group (Figure 8a). However, the more suitable area for the distribution of the Japan group was only in the southern and central Japan (Figure 9a). The northern boundary for the potential population distribution area of the Japan group in the East China Sea was Zhoushan. Climate change is predicted to influence the two groups differently. Future prediction for the China group under scenario RCP 4.5 showed that habitat suitability would clearly increase from south to north and the potential distribution may shift northward in the coastal waters of China and Japan by 2100 (Figure 8b,8c). However, the predicted result revealed a reduction in the potential area for the Japan group and no clear increase in the northern Japan areas by 2100. (Figure 9b,9c) This potential reduction in habitat suitability may indicate lower evolutionary potential to adapt to a changing environment. The response curves to temperature suggested different temperature adaptations between the two groups. The China group prefers high temperature, whereas the Japan group favors low temperature (Figure 8d,9d).