Genome‐wide evolutionary signatures of climate adaptation in spotted sea bass inhabiting different latitudinal regions

Abstract Consideration of the thermal adaptation of species is essential in both evolutionary biology and climate‐change biology because it frequently leads to latitudinal gradients of various phenotypes among populations. The spotted sea bass (Lateolabrax maculatus) has a broad latitudinal distribution range along the marginal seas of the Northwest Pacific and thus provides an excellent teleost model for population genetic and climate adaptation studies. We generated over 8.57 million SNP loci using whole‐genome resequencing from 100 samples collected at 14 geographic sites (five or ten samples per site). We estimated the genetic structure of the sampled fish and clustered them into three highly differentiated populations. The genetic differentiation pattern estimated by multivariable models combining geographic distance and sea surface temperature differences suggests that isolation by distance and isolation by environment both have significant effects on this species. Further investigation of genome‐wide evolutionary signatures of climate adaptation identified many genes related to growth, muscle contraction, and vision that are under positive natural selection. Moreover, the contrasting patterns of natural selection in high‐latitude and low‐latitude populations prompted different strategies of trade‐offs between growth rate and other traits that may play an essential role in adaptation to different local climates. Our results offer an opportunity to better understand the genetic basis of the phenotypic variation in eurythermal fishes inhabiting different climatic regions.


| INTRODUC TI ON
Water temperature is one of the most critical climatic factors for fishes, determining various physiologic phenotypes and strongly affecting their fitness and distribution (Payne et al., 2016). The natural population structure of fishes is also sensitive to water temperature due to adaptation to local climates (Sandstrom et al., 1995). Many studies have suggested that fishes living in different temperature environments undergo spatially heterogeneous natural selection (Bradbury et al., 2010;Hasselman et al., 2013;Lannig et al., 2003).
Many ecological studies have revealed that teleosts that are distributed at high latitudes tend to have higher growth potential (i.e., a higher maximum growth rate) than their low-latitude conspecies (Conover et al., 1997;Imsland et al., 2000;Lindgren & Laurila, 2005).
High growth potential is usually accompanied by reductions in resistance to predation , swimming performance , and metabolic rate (Sylvestre et al., 2007), which can be described as different trade-off strategies in resource allocation (Steiner & Pfeiffer, 2007). Indeed, adaptation to different temperatures involves shifts in a broad range of growth-related trade-offs. For example, Roze et al. (Roze et al., 2013) identified tightly interrelated trade-offs among thermal growth sensitivity, hypoxia tolerance, and growth in rainbow trout (Oncorhynchus mykiss). Matte et al. (Matte et al., 2020) studied the mortality-to-growth ratio of neighbouring brook trout (Salvelinus fontinalis) populations and revealed that this ratio increased linearly with increasing temperature. This could result in different trade-offs related to mortality and growth in this species. In addition, it has been widely discussed in recent years that trade-offs between growth and predator avoidance involving various physiological and behavioural traits such as swimming ability (Watson et al., 2018), lipid allocation (Mogensen & Post, 2012) and surplus aerobic metabolic capacity (Norin & Clark, 2017) play a critical role in teleost climate adaptation (Mogensen & Post, 2012;Norin & Clark, 2017;Watson et al., 2018).
The spotted sea bass (Lateolabrax maculatus) was recently discriminated from the Japanese sea bass (L. japonicus) as a newly recharacterized species with noticeable morphological and genetic differences (Yokogawa, 1998). Taking advantage of their eurythermic and euryhaline characteristics, spotted sea bass are distributed from the tropical Beibu Gulf in the South China Sea to the midtemperate Bohai Gulf (Figure 1a). The sea surface temperature (SST) in winter can fall to below 2°C in the northern Bohai Gulf and remain as high as 15°C in the Beibu Gulf, while the difference in SST in these two locations is very small in summer. Spotted sea bass inhabiting high-and low-latitude areas show different life-history strategies, as do many ectothermic animals. The northern populations in the Bohai Gulf have stronger growth potential and motor ability and higher cold resistance than their southern conspecifics in the South China Sea (Hu et al., 2019), making this species a good model through which to define the genetic basis of thermal adaptation. Our previous population genetic study employing double-digest restriction site-associated DNA (ddRAD) technology preliminarily investigated the population structure of this species (Zhao et al., 2018).
However, the lack of a reference genome made it impossible to obtain a sufficient number of SNP loci to support an in-depth investigation of the genetic mechanism that underlies thermal adaptive evolution in this species.
In the current work, we used a whole-genome sequencing (WGS) approach to detect genome-wide SNP sites based on an assembled reference genome of L. maculates (Chen et al., 2019) and reinvestigated the population structure in this species. Furthermore, we investigated how geographic and climatic factors drove the genetic differentiation of this species using isolation-by-distance (IBD) and isolation-by-environment (IBE) approaches and used a combination

| Sampling, sequencing, and SNP calling
We performed whole-genome resequencing of 100 spotted sea bass collected from 14 sites (ten samples from each site in the Beibu Gulf and the Bohai Gulf and five samples from each other site) along the coastline of China ( Figure 1a and Table S1). White muscles were collected from fish anesthetized with eugenol (200 mg/L for 15 min).
Muscle samples were lysed in SDS digestion buffer containing proteinase K. DNA extraction was performed using the phenolchloroform method. For each sample, a 250-bp insert-size library was prepared from approximately 0.5 μg of DNA according to the procedures outlined in the Illumina DNA Prep Reference Guide. All libraries were sequenced on the Illumina HiSeq 4000 platform, generating 150-bp paired-end reads.
Our previous reference genome assembly was completed with high accuracy and contiguity at the scaffold level (GenBank accession: GCA_004028665.1) (Chen et al., 2019). This version of the genome contains 1639 scaffolds ranging in size from 952 bp to 12.09 Mbp and has an N50 size of 2.79 Mb. However, LACHESIS software, which is now obsolete, was used to process the raw Hi-C data and chromosome-level scaffolding, and this may have introduced suboptimal orderings and orientations of the scaffolds . Hence, the genome was first refined using a better Hi-C data processing tool, 3D-DNA (version 201008) (Dudchenko et al., 2017). The improved genome assembly, which contained 480 scaffolds and had a total length of 589.99 Mbp, had higher integrity than the previous assembly ( Figure S1 and Table S1 in Appendix S1).
The improved genome also shows higher collinearity with the linkage map of L. maculatus and the genomes of related species (Figures S2 and S3,and Table S2 in Appendix S1). The quality of the raw sequencing data was assessed and controlled using SolexaQA++ (version 3.1.1) (Cox et al., 2010). SNP calling was then performed using GATK (version 4.0.2.1) (McKenna et al., 2010) under the guidance of GATK Best Practices (van der Auwera et al., 2013). Finally, hard filtering was applied to the raw SNP callset. Detailed information about these steps can be found in Appendix S1.

| Population genetics analysis
As highly linked SNPs are redundant for population genetics analysis, the final set was shrunk using PLINK with the parameter "--indep-pairwise 5 1 0.9" in each population. Three genetic diversity indices, observed heterozygosity (Ho), gene diversity (Hs), and Fis statistics, were calculated using the R package "hierfstat" (version 0.5-7) (Goudet, 2005). The mean Tajima's D across all chromosomes was estimated using "VCFtools" (version 0.1.17) (Danecek et al., 2011). In principal component analysis (PCA), the original variant call format (VCF) file was converted to CoreArray Genomic Data Structure (GDS) format using the R package "gdsfmt" (version 1.22.0) (Zheng et al., 2012). Eigenvectors and eigenvalues were then calculated using the R package "SNPRelate" (version 1.16.0) (Zheng et al., 2012) with default parameters. "Frappe" software (version 1.0) (Tang et al., 2005) (version 1.1) was employed to estimate the genetic ancestry of each sample with 100,000 expectation maximization (EM) iterations and up to four ancestral populations.
We constructed a neighbour-joining (NJ) tree of all samples using a command-line-based collection of utilities, VCF-kit (version 0.1.6) (Cook & Andersen, 2017). The size histories of all populations were estimated from whole-genome sequence data by SMC++ (version 1.15.4) (Terhorst et al., 2017). As molecular clocks have not been calculated in L. maculatus or its closely related species, a mutation rate of 2.0 × 10-9 per base per generation calculated in the Atlantic herring (Clupea harengus) (Feng et al., 2017) was used in this study.
Homozygous regions longer than 5 kb were treated as missing.
The polarization error was set to the default value (0.5). All other arguments were set to default values.

| Isolation by distance and environment analysis
The pairwise F ST values among populations were calculated within each chromosome using Plink2 (version v2.00a3) (Chang et al., 2015) with the "--fst" parameter based on Weir and Cockerham's method (Weir & Cockerham, 1984). Genetic distance was calculated as F ST / (1 − F ST ) (Rousset, 1997). Three geographic distances (great-circle distance, D gcc ; distance along the coastline, D csl ; and distance along latitude, D lat ) and three temperature differences (annual maximum, minimum, and median sea surface temperature differences, denoted SST max , SST min , and SST min , respectively) were used to model the genetic differentiation. D gcc was calculated from longitude and latitude coordinates using the following equations: where R, φ1, φ2, and ∆λ are the Earth's radius (6371.393 km), the latitude of one site, the latitude of the other site, and the difference in longitude between the two sites, respectively. D csl was calculated based on the preprocessed OpenStreetMap data obtained from the Planet OSM website (https://osmda ta.opens treet map.de). Raw data in shapefile format were transformed into "Mercator" projection using the R package "sf" (version 0.9.7) and then divided into 0.1 × 0.1-km 2 grids. Finally, D csl was calculated as the sum length of the diagonal lines of all grids connecting two places. Daily SSTs from 1971 to 2010 were extracted from NOAA high-resolution SST data (https://psl.noaa.gov/) (Reynolds et al., 2007). The "netCDF4" Python library (version 1.5.5.1) was used to parse the NetCDF4 file. We averaged SST values within a 0.25-degree grid covering each site.
To model how a single factor influences the genetic distance, single-factor linear regression between genetic distances and exogenous geographic and climatic variables was implemented using the "statsmodels" Python library (version 0.12.1) for IBD and IBE analyses. When attempting to dissect the effects of IBD and IBE patterns, the main challenge we faced was multicollinearity among exogenous variables. To address this, three mutually independent methods with high tolerances for multicollinearity were employed to solve the problem. Regularized least squares (RLS) regression is a family of relatively new linear regression techniques that were recently introduced for solving problems of confounding and causality plaguing genome-environment association studies (Frichot & Francois, 2015). The generalized dissimilarity model (GDM) is a matrix regression technique that is powerful when applied to the model genetic distance between sampling sites against geographic and environmental variables (DeSilva & Dodd, 2020; Geue et al., 2016).
The third method, canonical correlation analysis (CCA), is a powerful and widely used method for coping with multicollinearity in IBD and IBE analysis (Prunier et al., 2015). We conducted regularized leastsquares regression using the "linear_model" module contained in the "sklearn" python library (version 0.24.0). In RLS regression, unimportant variables were first discarded by least absolute shrinkage and selection operator (LASSO) regression. The remaining variables were then used to build an elastic net model with an alpha value of 0.5. Both LASSO and elastic net models were built using the "sklearn" python library (version 0.24.0). The "gdm" R package (version 1.4.0) (Mokany et al., 2022) was employed to fit a GDM. The model significance and variable importance were estimated using the "varImp" function with 5000 matrix permutations. The "isplineExtract" function in the "gdm" R package was used to extract an I-spline whose maximum height indicates the significance of predictors (Ferrier et al., 2007) for each predictor from the resulting GDM objects. In CCA, all descriptive statistics were calculated using the "calc.yhat" function in the "yhat" R package (version 2.0-3) (https://cran.uib.no/ web/packa ges/yhat/index.html). Their confidence intervals were generated after 5000 bootstraps by the "booteval.yhat" function.

| Identification of temperature-associated SNP loci and genes
Temperature-associated SNP loci were identified based on SST max and SST min using Baypass (version 2.2) (Gautier, 2015). The population covariance matrix was first estimated in the core model with default parameters and was then used as an extra input in standard covariate mode to estimate the Bayes factor (BF) of each SNP with default parameters. Bayes factors were estimated in deciban (dB) units. Only SNPs with a decisive evidence level (BF > 20 dB) of associations were identified as SST max -or SST min -associated SNP loci (Hoijtink et al., 2019). RDA was performed using the "vegan" (version 2.6-4) library in R (Dixon, 2003). The RDA model was first built using the "rda" function with default parameters. The p values of SNPs were then calculated from their loading scores with a chi-square distribution. The overlaps between outliers identified using BayPass and vegan were considered temperature-associated SNPs. Finally, Metascape (version 3.5)  was employed for enrichment analysis based on genes located around these SNPs. This software was also employed to calculate nucleotide diversity values. The ratio of nucleotide diversity in the BH population to that in the BB population was calculated and then log-2 transformed

| Selective signature detection
(log 2 π BH /π BB , denoted ω). The extent of haplotype homozygosity (Rsb) was calculated for each SNP using the R package "rehh" (version 3.0.1) (Gautier et al., 2017) and then averaged within sliding windows. Positively selected regions (PSRs) were identified as regions having absolute F ST , ω and Rsb in the 95th percentiles. Genes whose gene bodies overlapped with any PSR were identified as positively selected genes (PSGs). In addition, Tajima's D index was estimated using VCF-kit (version 0.1.6) (Cook & Andersen, 2017) with the parameter "tajima". The composite selection score (CSS) (Avalos et al., 2017)

| Prediction of the tertiary structure and thermostability of myosin heavy chain (myhz2)
To understand the functional impact of the identified mutations in the myosin heavy chain (myhz2) gene on protein thermostability, the dominant alleles of each population and each nonsynonymous mutation were first applied to the reference myhz2 coding sequence to derive BH-and BB-specific protein sequences. The tertiary structures assumed by proteins with these two sequences were then predicted using I-TASSER (version 5.1) (Yang et al., 2015) with an identity cut-off of 0.3. Two approaches were employed to estimate the thermostability of the two structures. The "SCooP" web server (version 1.0) (Pucci et al., 2017) was used to predict the Gibbs-Helmholtz equation associated with the folding transition. In addition, another web server, deepDDG (Cao et al., 2019), was used to predict the changes in protein stability brought about by point mutations using neural networks.

| Sequencing and SNP detection
High-throughput sequencing generated 3.97 billion pairs of raw reads with a total length of 1.19 Tbp and an average depth of 19.64×.
After quality control, 98.42% of the reads (1.17 Tbp) remained for processing in the downstream steps. Of these, 783.6 billion clean reads were mapped to the reference genome of L. maculatus at mapping ratios ranging from 97.2% to 99.1%. We obtained a total of 8.57 million SNPs from the 100 sampled fish (Figure 1a, Figure S1, and Tables S1-S3). We also generated a pruned SNP set after the removal of 1.19 million (13.88%) highly linked SNPs. This pruned set was used only in the genetic structure analysis.

| Genetic structure of the spotted sea bass
As shown in Figure 1b implying that they are highly differentiated. Phylogenetic trees support the same division of these samples ( Figure 1c and Figure S2).
Ancestry proportion estimation reconfirmed that the samples can be divided into three populations when K is set to three. When K is set to two, the IM populations are depicted as an admixture of the BH and BB populations. Most of the samples in the IM population, with the exception of the ZJ samples, were found to harbor more genetic material introgressed from the BH population than from the BB population, (Figure 1d).
We built an NJ tree by taking a related species (the Asian sea bass, AS) as the outgroup (Figure 1c). The rooted phylogenetic tree showed that the spotted sea bass found in Chinese coastal seas originated in the Bohai Gulf and then spread to southern areas.
In addition, the historical effective population sizes suggest that during this process the BB populations experienced rapid expansion 410,000 years ago followed by a bottleneck occurring approximately 290,000 years ago (Figure 1e).

| Isolation by distance and environment analysis
In the IBD analysis, models fitted using D csl and D gcc had higher goodness of fit (R 2 adj = 0.652 and 0.607, AIC = −540.943 and −531.416, for D csl and D gcc , respectively) and greater significance (p = 2.57E-19 and 2.76E-17) than those fitted using D lat (R 2 adj = 0.444, AIC = −504.385, and p = 1.67E-11) ( Figure 2a, Table 1, and Figure S3). In IBE analysis, SST max showed the sole negative Akaike information criterion (−32.899) and a p value (7.53E-16) that was much lower than those associated with the other two factors, indicating a more significant correlation with genetic distance (Figure 2b, Table 1, and Figure S4).
We also used three multivariable methods to dissect the effects of IBD and IBE patterns on spatial genetic divergence. This analysis was greatly complicated by the multicollinearity among geographic factors and climatic factors. Using these multivariable methods, we found that D csl and SST max simultaneously explained most of the variation in genetic distance and that SST max had a greater effect than did D csl on genetic distance. In LASSO regression, only D csl and SST max had nonzero coefficients (0.679 and 0.283, respectively, Table 2). Similarly, D csl and SST max were the only two factors retained in the final GDM (for these two factors, the residual deviance of the response variable was 3.941 and 12.469, respectively) ( Table S4). All indices produced by CCA indicated that D csl and SST max were the two most significant predictors ( Figure 2c and Table S5). Although the Pratt measure and beta weight suggested that D csl had higher importance than SST max , all other predictors revealed a more significant impact of SST max (Figure 2d and Table S5).
We then divided the genome into PSRs and non-PSRs (see Sections 2.6 and 3.5) and rebuilt a GDM for each division. The non-PSR GDM was very similar to the whole-genome GDM (Figure 2e,f).
However, the magnitude of the total genetic distance change along D csl and SST max showed very different patterns in PSR-specific GDM.

| Identification of temperature-associated SNP loci and genes
We combined a Bayesian approach (BayPass) with RDA to detect SNPs associated with SST max and SST min . BayPass revealed that 516 and 227 SNPs had strong associations solely with SST max or SST min , respectively, and that 164 SNPs were associated with both SST max and SST min (Figure 3a and Table S6). The RDA detected 1927 outlier SNPs. A total of 334 SNPs were detected with both RDA and BayPass and identified as candidate temperature-associated SNPs (178 SST max -associated SNPs, 20 SST min -associated SNPs, and 146 SNPs associated with both SST max and SST min ; Figure S5). We then searched genes located within 5 kbp of the candidate SNPs. In that search, we identified 186 SST max -associated genes and 79 SST minassociated genes (Figure 3a and Table S6). Functional enrichment revealed that the SST max -associated genes were significantly enriched in 6 functional clusters, the top two of which were closely related to nervous system function ( Figure S6 and Table S7). SST min -associated genes were enriched in only three functional clusters ( Figure S6 and Table S8). The enriched GO biological process, "muscle contraction" (GO:0006936), comprised four genes that encode components of muscle fibres and their binding proteins, proteins that have direct effects on locomotion. The observed enrichment of the "insulin signalling" pathway (WP1313) indicates that regulation of growth plays a role in the cold adaptation of BH populations

| Candidate genomic regions under extreme positive selection
While Bayesian approaches are powerful for identifying SNP loci that harbour subtle shifts in allele frequency associated with environmental gradients, such approaches have difficulty detecting different haplotypes associated with individual genes (Salmón et al., 2021). Genome-wide scanning for selective sweeps is a common and powerful way to identify candidate regions and genes that are subject to local adaptation and are involved in responses to climate change (Chen & Narum, 2021;Han et al., 2021;Li et al., 2021).
Using this method, 436 and 416 PSRs were identified in the BH and BB populations, respectively (Figure 3, Figure S7, and Table S9). Note: RLS models were built using LASSO regression and elastic net regression. The columns in this table show the geographic/climatic factors used to build the models and the coefficient, standard error, t statistic, and p value of each factor. Each model is also summarized according to the coefficient of determination (R 2 ), SSE, AIC, and BIC.
TA B L E 2 Regularized least squares (RLS) models indicating the relative importance of geographical and climatic factors.
Approximately one-third of the temperature-associated SNPs were located on these PSRs. However, none of them were located within 1000 bp of the nearest PSR. The other two-thirds of the temperature-associated SNPs were distributed away from PSRs at distances ranging from 1000 bp to 10 Mbp ( Figure S8). We noted that BB-PSRs appeared to have experienced more intense natural selection than did BH-PSRs. Stronger selective sweeps in BB-PSRs were indicated by ω, Rsb, Tajima's D, and CSS, whereas higher F ST suggested that the BH-PSRs are more thoroughly differentiated ( Figure 3e and Table S10). One reasonable explanation for this is that BB populations have a stronger but shorter selection history.
We identified 287 and 227 protein-coding genes as BH-and BB-PSGs, respectively (Tables S11 and S12). Within these genes, we identified 406 highly differentiated nonsynonymous SNPs (Table S13). Notably, many genes that encode essential structural components or regulatory proteins of striated muscle fibres were identified as PSGs in the BH population; these included myhz2, titin, and cofilin-2 (Figure 4a and Figure S9). We identified two nonsynonymous mutations (c.373C > A and c.700G > A) in myhz2 that were highly differentiated between the BH and BB populations ( Figure 4b). The former (c.373C > A) also has a significant association with SST min (Bayes factor = 21.14; Figure 3a and Table S6). Although the tertiary structure of the myosin heavy chain did not change significantly in the presence of this mutation (Figure 4c), the thermal stability of this protein decreased in the BH population ( Figure 4d and Table S14).
We conducted overrepresentation analysis (ORA) and subsequent enrichment clustering to explore biological processes  Table S15). Thirty PSGs in the southernmost BB populations were enriched in eight functional categories, and these were further classified into eight enrichment clusters.
Two relatively large clusters represented by the GO terms "ceramide biosynthetic process" (GO:0046513) and "embryonic eye morphogenesis" (GO:0048048) were directly related to the development and function of the visual system ( Figure 5d); the genes in these clusters included fgf8, vax1, pard3ab, foxg1, and ipo13. The two most significant clusters, which were represented by the GO terms "fibroblast growth factor receptor signalling pathway" (GO:0008543) and "embryonic neurocranium morphogenesis" (GO: 0048702), contained a series of genes that are also probably instrumental in visual system development (Figure 5e,f, and Table S16).

| Genetic structure and evolutionary history of spotted sea bass populations
The whole-genome sequencing conducted in this study generated a sampled spotted sea bass (BH, IM, and BB) that differ greatly in their genetic compositions (Figure 1). Moreover, we observed that the habitats of these three populations are separated by the Qiongzhou Strait rather than by the brackish water present in large rivers.
L. maculatus lives in nearshore waters and has limited migration ability (Sun et al., 2014) and wide salinity tolerance (Shao et al., 2018;Zhang et al., 2017); thus, it is less able to traverse geomorphological barriers such as peninsulas and shallow waterways than to spread across hydrochemical obstacles.
The Fish are less active in cold waters due to reduced muscle sensitivity caused by low temperatures (Riisgard & Larsen, 2009;Rome et al., 2000;Yan et al., 2012). In accordance with previous findings, many PSGs that encode essential structural components or regulatory proteins of striated muscle fibres were identified in the BH population ( Figure 4a and Figure S9), indicating that changes in swimming ability may be a key behavioural adaptation process in the spotted sea bass. Among these genes, myhz2 was repeatedly highlighted by selective signatures, association with SST min , and the presence of highly differentiated nonsynonymous mutations (Figures 3a and 4). Computational thermochemical analysis revealed that these mutations decrease the thermostability of the protein encoded by this gene. A previous study indicated that reduced thermostability of myosin is responsible for its more flexible structure in cold-acclimated carp and that this reduced thermostability further results in a reduced activation enthalpy for contraction at low temperatures (Watabe, 2002). Such a reduction in the thermal stability of myofibrillar proteins has been reported as an adaptive mechanism that allows teleosts and other vertebrates to cope with environmental temperature changes (Howell et al., 1991;Rodgers et al., 1987;Takahashi et al., 2005). Overall, the evolutionary signatures of sarcomere-related genes in BH populations of spotted sea bass revealed the presence of an adaptive mechanism for thermal acclimation that involves changes in intrinsic muscle contractile properties in teleosts.

| Northern and southern populations may make different trade-offs between growth and mortality
Our IGFBPs play a crucial role in the growth of the skeletal, muscular, and nervous systems in teleosts (Shimizu & Dickhoff, 2017). In addition, we identified selective signatures of the suppressor of the cytokine signalling-2 (socs-2) gene, an important negative regulator of the GH/IGF-1 axis. The product of this gene can competitively bind to Igf-1r and thereby impede the function of IGF-1 in various species (Greenhalgh et al., 2002;Horvat & Medrano, 2001;Liu et al., 2006;Philip & Vijayan, 2015).
In the BB population, both the Bayesian and selective signature scanning approaches identified PSGs that are closely related to the development and function of the nervous system (Figures 3a and 5d-f; Tables S7 and S16). Notable among them were genes related to visual sensing and functional categories (Figures 5d-f and S11; Table S16). Although other types of sensation are involved, visual perception plays a crucial role by allowing organisms to detect predators and prey (McCormick & Manassa, 2008;Oliveira et al., 2017;Puvanendran et al., 2002). In tropical zones, predator-prey interactions are stronger and more complex than those in colder areas (Schemske et al., 2009). In response to this, fishes generate a variety of alterations in their behaviour, physiology, and demography (Culumber, 2020;Mennen & Laskowski, 2018;Mitchell & Harborne, 2020;Xu et al., 2019).
The divergent patterns of PSG enrichment found in the two populations suggest that the populations may adopt different strategies that involve different trade-offs between growth and mortality. It has been demonstrated that a higher capacity for growth is a common adaptation to shorter growing seasons in high-latitude environments (Angilletta et al., 2002;Conover & Schultz, 1995;Suzuki et al., 2010). However, rapid growth is realized through frequent feeding and is usually accompanied by a high risk of predation, a situation that has been termed the "growth-mortality trade-off" Stocks et al., 2005;Suzuki et al., 2010). A common garden experiment demonstrated that northern populations of L. maculatus had a significantly higher growth rate (2.4 g/d) than their southern conspecifics (1.7 g/d) (Chen et al., 2001). However, the presence of a latitudinal difference in predator avoidance remains to be validated in this species.
It is worth noting that dissection of the interactions of the genome of an organism with the environment is very complicated. The presence of relatively sparse sampling sites and multicollinearity among geographic and environmental factors undoubtedly increase the challenge. Therefore, the findings of this study should be validated by quantitative studies of adaptive traits followed by common garden-variety experiments in the future.

ACK N OWLED G M ENTS
We acknowledge financial support from the Guangdong Basic

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare that they have no conflict of interest. The authors alone are responsible for the content and writing of the article.