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).